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I.  INTRODUCTION 


A.  BACKGROUND 

1.  A  Carrier  Navy 

The  carrier  is  the  nucleus  of  today's  U.S.  Navy,  but  its  ability  to  project  power 
is  a  byproduct  of  the  synergistic  relationship  it  shares  with  its  air  wing.  As  Peter 
Garrison  stated  in  CV:  Carrier  Aviation,  "by  adopting  submarine-launched  ballistic 
missiles  and  ship-based  tactical  aircraft  the  navy  has  co-opted  some  of  the  tools  of  the 
army  and  air  force  and  adapted  them  to  a  marine  environment.  Ships  in  themselves  are 
no  longer  powerful,  they  take  their  power  from  the  airplanes  or  long  range  missiles  that 
they  carry.  An  adversary  setting  out  to  sink  a  ship  does  not  send  another  ship  against 
it,  but  instead  sends  an  airplane  or  a  missile"  [Ref.  1].  (italics  mine).  The  fleet  has 
become  dependent  upon  aircraft  to  take  the  battle  over  the  horizon,  and  to  extend  the 
battle  group's  area  of  influence.  Without  downplaying  the  importance  of  the  carrier,  it 
is  just  another  large  ship  without  its  contingent  of  aircraft.  Thus,  the  carrier's  real 
importance  to  the  feet  rests  on  the  ability  of  its  personnel  to  support  and  maintain  the 
air  wing  at  an  optimal  level  of  readiness. 

2.  Aircraft  Readiness 

Ideally,  the  fleet  would  like  each  aircraft  to  be  fully  mission  capable  (FMC); 
that  is,  the  material  condition  of  the  aircraft  should  satisfy  minimum  requirements  to 
perform  all  of  its  missions.  An  aircraft's  mission  status  (e.g.,  partially  mission  capable 
(PMC),  not  mission  capable  (NMC),  etc.)  is  determined  using  the  Mission  Essential 
Subsystems  Matrices  published  in  OPNAVINST  5442. 4H. 

The  ability  to  ensure  a  high  percentage  of  FMC  aircraft  is  constrained  by 
several  factors.  First,  the  Navy  has  not  been  able  to  provide  shipboard  intermediate 
maintenance  activities  (IMA)  with  sufficient  numbers  of  qualified  personnel,  especially 
in  critical  technical  billets.  In  part,  this  results  from  the  allocation  of  manpower  based 
on  squadron  size  (i.e.,  the  smaller  squadrons  have  fewer  personnel  to  maintain  their 
aircraft).  The  problem  associated  with  the  manpower  shortage  is  compounded  by  a 
need  for  specializatiop  in  troubleshooting  and  repairing  the  different  aii craft  types, 
models,  and  ^scries  deployed  on  the  carrier.  Second,  the  IMA's  capability  to  repair 
components  as  they  fail  is  constrained  by  limitations  of  test  equipment,  maintenance 
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facility  size,  availability  of  repair  parts,  and  the  total  demand  within  the  repair  pipeline. 
Finally,  those  components  which  cannot  be  fixed  at  the  intermediate  level  must  be 
resupplied  by  outside  sources. 

The  air  wing  must  be  able  to  accurately  identify  and  communicate  its 
requirements  for  component  replacements  and  resupply  facilities  should  be  capable  of 
expeditiously  responding  to  the  those  needs.  During  periods  of  high  tempo  flight 
operations,  as  would  be  experienced  during  wartime,  the  factors  listed  above  become 
more  critical.  The  ability  to  maintain  high  levels  of  readiness  is  hampered  by  an 
increase  in  the  number  of  components  in  the  repair  pipeline,  manpower  and  facility 
limitations,  and  a  reduced  or  nonexistent  ability  to  resupply  vital  components. 
Potential  problems  may  be  partially  rectified  by  improving  the  methodologies  currently 
used  to  attain  that  "best"  possible  mix  of  replacement  parts  located  cn  the  aircraft 
carrier.  [Ref.  2:  pp.  3,9] 

B.  AVIATION  CONSOLIDATED  ALLOWANCES 

1.  Definition 

The  Aviation  Consolidated  Allowance  List  (AVCAL)  is  supposed  to  be  a 
"best"  mix  inventory.  An  AVCAL  is  defined  in  OPNAVINST  479G.2C  as  a  list  of 
aeronautical  material  tailored  to  each  individual  ship  to  support  assigned  or  embarked 
aircraft  flight  operations.  This  list  of  parts  is  prepared  by  the  Navy's  Aviation  Supply 
Office  (ASO)  under  the  direction  of  the  aircraft  type  commander,  and  determines  spare 
part  allowances  based  on  expected  demands. 

Ideally,  it  would  be  advantageous  to  stock  at  least  one  spare  for  each  part 
that  could  possibly  fail  during  deployment  regardless  of  how  unlikely  that  occurrence 
might  be.  Unfortunately,  inventory  sizes  are  constrained  by  both  budget 
considerations  and  storage  limitations.  The  latter  is  the  result  of  a  diverse  collection  of 
aircraft  types,  manufactured  by  different  companies  with  each  aircraft  requiring  its  own 
specific  replacement  parts,  coupled  with  the  inherent  aversion  to  building  aircraft 
carriers  even  larger  than  they  already  arc  simply  to  carry  more  parts.  Table  1  lists  the 
normal  complement  of  aircraft  aboard  a  carrier,  and  is  provided  to  illustrate  the 
differences  in  aircraft  types,  missions,  and  spare  part  requirements.  [Ref.  2:  p.5). 

2.  Need  for  Improvement 

Determination  of  the  proper  AVCAL  allowance  for  each  component  is  not  an 
easy  task,  and  is  definitely  an  area  needing  improvement.  A  recent  study  conducted  by 
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TABLE  1 

NORMAL  COMPLEMENT  OF  CARRIER  AIRCRAFT 


AIRCRAFT 

MISSION 

SQUADRONS 

AIRCRAFT 

F-14A 

Fighter 

2 

24 

A-7E 

Light  Attack 

2 

24 

A-6E 

Attack 

1 

10 

KA-6D 

Tanker 

1 

4 

S-3A 

ASW 

1 

10 

E-2C 

SURV 

1 

4 

EA-6B 

EW 

1 

4 

SH-3 

ASW 

1 

6 

TOTAL 

81 

the  Center  for  Naval  Analyses  reported  that  statistics  for  combat  aircraft  operating 
during  the  period  July  1982  through  June  1983  show  aircraft  operating  at  less  than  60 
percent  of  their  anticipated  wartime  rate.  These  aircraft  were  not  FMC  due  to  supply 
between  18-24  percent  of  the  time.  These  readiness  rates  would  drop  significantly  if  the 
aircraft  flew  at  a  wartime  rate  and  AVCALs  continued  to  be  constructed  as  they  were 
in  this  period"  [Ref.  3:  p.1-1].  This  statement  strengthens  the  argument  that  current 
methods  of  AVCAL  formulation  are,  at  best,  adequate  for  peacetime  operations,  but 
could  result  in  disastrously  low  readiness  rates  wnen  exposed  to  the  additional 
requirements  associated  with  wartime  utilization.  Two  approaches  to  determining  a 
better  AVCAL  would  be  the  introduction  of  new  methods  of  modelling  requirements  or 
the  improvement  of  weak  areas  within  the  current  model. 

3.  AVCAL  Models 

Two  categories  of  models  used  to  compute  the  AVCAL  are  readiness-based 
and  demand-based  inventory  systems.  Readiness-based  systems  determine  spare 
allowances  that  achieve  a  predetermined  level  of  readiness  while  minimizing  the  cost  of 
the  AVCAL.  Examples  of  readincss-bascd  models  arc  Rand's  DYNAMETRIC  model 
and  the  Center  for  Naval  Analyses  MIME  model  [Ref.  3:  P-2-8J.  Demand-based 
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(supply  effectiveness)  models  determine  component  allowances  based  on  expected 
demand  and  the  capabilities  of  the  repair  facility  and  resupply  pipeline.  The  ASO 
Manual  model  is  an  example  of  a  demand-based  method,  and  is  the  model  currently 
used  by  the  Navy  to  compute  AVCAL  spare  allowances. 

C.  ASO  MANUAL  MODEL 

The  ASO  Manual  Model  uses  historic  data  to  forecast  demand  rates  and 
component  repair  times.  The  model  determines  the  type  and  number  of  components 
included  in  the  AVCAL  allowance.  For  a  detailed  explanation  of  the  rules  governing 
component  selection,  the  reader  should  refer  to  either  Aviation  Parts  Allowance  Policy 
by  Peter  Evanovich  or  A  Retail  Inventory  Model  for  Saval  Aviation  Repairable  Items  by 
Mark  L.  Mitchell.  [Refs.  3,4]. 

1.  Attrition  Only  Items 

Attrition  only  items  are  repairable  components  that  cannot  be  repaired  at  the 
intermediate  maintenance  activity  aboard  the  carrier.  The  expected  demand  for  items 
of  this  type  is  computed  using  equation  1.1. 

NB'-NB  x(FH'.'FH)  (eqn  l.l) 

where, 

NB'  “  number  of  expected  attritions  for  a  specific  flight  program 

NB  =  number  of  attritions  in  the  ASO  data  base 

FH'  =  specified  flight  program 

FH  =  total  (light  hours  represented  by  the  ASO  data  base 

2.  Items  that  can  be  Repaired  Locally 

The  allowance  for  items  that  can  be  repaired  locally  is  composed  of  two  parts. 
An  allowance  for  repairable  components  that  are  beyond  the  capabilities  of  the  local 
maintenance  activity  to  repair  them  (BCM),  and  an  allowance  for  components  that  are 
queued  at  various  places  in  the  local  repair  cycle. 

a.  BCM  Allowance 

Units  are  considered  BCM  if  the  maintenance  activity  cannot  repair  the 
item  because  of  a  lack  of  tools,  equipment,  parts  or  technical  expertise.  Equation  1.1  is 
used  to  forecast  the  number  of  BCM  failures  expected  to  occur  for  a  specific  flight 
program. 
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b.  Repair  Pipeline  Allowance 

The  second  input  to  the  local  repairable  allowance  provides  units  to 
compensate  for  the  number  of  failed  components  undergoing  repair  at  the  IMA.  This 
allowance,  called  the  Local  Repair  Cycle  Allowance  (LRCA),  is  set  at  a  quantity  which 
should  ensure  that  demand  does  not  exceed  available  supply  at  least  ninety  percent  of 
the  time.  This  safety  level  is  dependent  upon  an  accurate  estimate  of  the  number  of 
components  repaired  by  the  shipboard  intermediate  maintenance  activity.  Such  a 
forecast  is  expressed  in  equation  1.2. 

NR'  =  NR  x(FHVFH)  (eqn  1.2) 


where, 

NR'  =  number  of  expected  repairs  for  a  specific  flight  program 
NR  28  number  of  repairs  represented  in  the  ASO  data  base 
”  and, 

FH'  and  FH  remain  as  defined  above 

D.  FAILURE  RATE  DISTRIBUTION 

Combining  equations  1.1  and  1.2  provides  the  expected  number  of  total  failures, 
both  repairables  and  attntables,  for  the  specified  flight  program.1 

NR'  +  NB'  -  (NR+NB)x(FH*;FH)  (eqn  1.4) 

Substituting  NF  for  NR'  +  NB'  and  NF  for  NR  +  NB,  equation  1.4  can  be  expressed 
as: 


NF'  *  NF  x  (FH  FH)  (eqn  1.5) 

Equation  1,6  results  from  rearranging  terms  on  the  right  hand  side  of  equation  1.5. 

NF'  =  (NF  FH)  x  FH'  (eqn  1.6) 


JFor  "attrition  only  items"  NR'  equals  zero. 
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The  total  number  of  expected  failures  NF‘  equals  the  product  of  some  specific 
future  flight  program  (FH')  and  the  quantity  NF/FH.  The  latter  should  be  easily 
recognized  as  the  maximum  likelihood  estimate  (based  on  the  historical  data  base)  for 
the  failure  rate  lambda.  Thus,  NF'  is  the  expected  value  of  a  Poisson  distribution  with 
the  failure  rate  parameter  X  and  interval  of  observation  (0,FH'].  The  use  of  the  Poisson 
distribution  implies  that  the  underlying  distribution  for  the  time  between  failures  is 
exponential.  Certain  properties  of  the  exponential  distribution  simplify  the 
computational  problem,  but  what  price  is  paid  for  this  simplicity?  The  remainder  of 
this  paper  will  suggest  alternative  distributions  for  describing  the  time  between 
component  failures  and  estimating  the  expected  number  of  failures  for  a  given  flight 
program. 

E.  SYNOPSIS 

Chapter  II  describes  the  methods  and  assumptions  used  to  extract  component 
failure  information  from  a  data  base  containing  six  squadrons  of  F-14A  aircraft  flight 
and  maintenance  records.  Chapters  III  and  IV  discuss  current  and  possible  methods 
of  modeling  aircraft  failure  rates,  and  techniques  used  to  derive  maximum  likelihood 
estimates  of  model  parameters.  Flight  hour  dependent  models  are  presented  in  Chapter 
Ill  and  include  the  traditional  Exponential  model  and  a  Weibull  modei.  Chapter  IV 
discusses  a  sortie-dependent  Geometric  model,  and  the  sortie-flight  hour  Mixed  model 
proposed  in  Distinguishing  the  Effects  on  Failures  of  Changes  in  Sortie  Rate  and  Sortie 
Length  by  Robert  A.  Levy  [Ref.  5],  The  four  ditferenl  models  will  be  compared  in 
Chapter  V  using  both  graphical  and  quantitative  measures  of  fit.  Finally,  Chapter  VI 
will  discuss  the  advantages  and  disadvantages  of  each  model,  along  with  the 
conclusions  and  recommendations. 


II.  DATA  BASE 


A.  ORIGIN  OF  DATA 

The  initial  data  base  was  obtained  from  the  Navy  Maintenance  Support  Office 
(NMSO),  Mechanicsburg,  Pennsylvania  via  the  Center  for  Naval  Analyses,  Alexandria, 
Virginia.  The  data  set  included  the  historical  flight  and  maintenance  records  for  two 
operational  deployments  of  F-14A  aircraft.  The  most  recent  deployment  included  two 
squadrons,  VF-114  and  VF-213,  which  operated  off  the  U5>S  ENTERPRISE  (CVN-65) 
during  the  period  of  March  1984  through  December  1984.  This  deployment  became  the 
"model"  set,  with  its  flight  and  maintenance  records  providing  the  parameter  estimates 
for  the  study. 

The  second,  earlier  deployment  consisted  of  four  squadrons,  VF-11  Y'F-31  and 
VF-14/VF-32,  which  operated  off  the  L'SS  KENNEDY  (CV-67)  and  I'SS 
INDEPENDENCE  (CV-62),  respectively,  during  the  period  of  November  1983  through 
March  1984.  In  thri  studv  these  four  squadrons  were  treated  as  the  "validation  '  set. 

Using  the  parameter  estimates  from  the  model  set,  four  specific  failure  rate 
distributions  were  used  to  calculate  the  expected  number  of  component  failures  which 
should  occur  during  a  deployment  similar  to  the  one  flown  by  the  four  squadrons 
aboard  the  USS  KENNEDY  and  L'SS  INDEPENDENCE.  The  actual  number  of 
component  failures  observed  during  the  deployment  of  these  squadrons  was  then 
compared  with  the  expected  number  generated  using  the  parameters  estimated  from  the 
"modei"  set  deployment. 

B.  MAINTENANCE  DATA  SYSTEM  (MDS) 

1.  Purpose 

The  maintenance  data  system  (MDS)  was  developed  to  improve  the 
management  capabilities  of  the  Navy's  Maintenance  and  Material  Management  (3-M) 
system.  It  was  designed  as  a  management  information  system  which  would  provide  the 
sta  .stical  data  needed  to  analyze: 

•  Equipment  maintainability  and  reliability 

•  Equipment  configuration  to  include  alterations  and  technical  directive 
compliance 

•  Equipment  mission  capability  and  utilization  rates 

•  Material  usage 
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•  Maintenance  and  material  procurement  times 

•  Weapon  system  and  material  costing 

To  meet  the  data  requirements  listed  above,  NMSO  divides  the  data  collection 
elTort  into  four  specific  areas:  material  reporting  (MR),  subsystem  capability  impact 
reporting  (SCIR),  maintenance  data  reporting  (MDR),  and  utilization  reporting.  The 
material  reporting  system  documents  those  materials  used  in  support  of  the 
maintenance  efforts.  The  subsystem  capability  impact  report  provides  the  higher  level 
commanders  with  information  pertaining  to  their  subordinate  command's  ability  to 
conduct  their  mission.  The  last  two  areas  of  data  collection,  maintenance  data 
reporting  and  utilization  reporting,  provided  the  data  necessary  to  determine  the  times 
between  component  failures. 

2.  Maintenance  Data  Reporting  (MDR) 

The  documents  which  support  the  maintenance  data  reporting  system  include 
support  action  forms  (SAF),  metrology  equipment  recall  (METER)  cards,  and  the 
Visual  Information  Display  System' Maintenance  Action  Forms  (VIDS  MAF).  Support 
action  forms  document  maintenance  actions  that  did  not  require  any  corrective  action, 
such  as.  aircraft  servicing,  aircraft  handling,  and  the  preliminary  look  phase  of  an 
inspection.  METER  cards  identify  the  maintenance  time  spent  calibrating  and  repairing 
the  test  and  measuring  systems.  The  last  form.  VIDS  MA  7  (OPNAV  4790  60), 
documents  inspections,  technical  directive  compliance,  and  repair  actions.  It  identifies 
the  component  failures  by  work,  unit  code,  time  of  failure  (julian  date  and  when 
discovered  code),  aircraft  experiencing  the  failure  (bureau  number),  cause  of  failure 
(malfunction  description  code),  and  type  of  maintenance  required  (transaction  code  and 
action  taken  code). 

3.  Utilization  Reporting 

The  flight  records  (utilization  rates)  were  compiled  from  information  recorded 
on  the  Aircraft  Yellow  Sheet  (OPNAV  3760, 2B).  The  Yellow  Sheet  is  prepared  by  air 
crew  personnel  at  the  end  of  each  mission  and  contains  information  identifying  the 
flight  according  to  aircraft  flown,  mission  type,  julian  date,  flight  duration,  etc. 

C.  SELECTED  WORK  UNIT  CODES 

A  work  unit  cede  (WLC)  is  an  alphanumeric  code  identifying  an  item  on  which 
work  is  being  performed.  Systems  are  identified  with  two  digit  codes  and  components 
with  five  digit  codes.  Seven  digit  codes  are  used  to  further  breakdown  components  into 
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lower  levels  of  subassembly.  The  components  selected  for  this  study  are  listed  in 
Table  2.  Each  of  the  nine  components  chosen  had  at  least  thirty  failures  during  the 
deployment  cycle. 


TABLE  2 

COMPONENTS  WITH  ASSOCIATED  WORK  UNIT  CODES 


COMPONENT  WUC 

CN 1 263/ ASN92( V)  INERTIAL  MEASURING  UNIT  734H100 

ECU74/A  SIGNAL  DATA  CONVERTER  46X1600 

MXU61 1; A  i:UEL  TANK  RELEASE  MECHANISM  4622100 

CPI  166  A  AIR  DATA  COMPUTER  56X2500 

AN  ARC159  159(V)-1  UHF  TRANSMITTER  632Z1XK) 

T1224  AWG9  RADAR  TRANSMITTER  74A1500 

IP11S5/AWG9  DETAIL  DATA  DISPLAY  74A5M00 

CP124S/ASW43  AIRCRAFT  ROLL  COMPUTER  5772200 

IP10Z7/AVA12  ANALOG  DISPLAY  INDICATOR  •  691S100 


I 
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D.  DATA  REDUCTION 

The  intent  of  the  data  reduction  process  is  to  eliminate  those  maintenance 
actions  which  do  not  identify  a  component  failure  (e.g.,  aircraft  servicing, 
troubleshooting,  calibrations,  or  other  routine  maintenance  actions).  The  maintenance 
codes  not  indicative  of  a  component  item  failure  are  found  within  the  series  of 
malfunction  description,  transaction,  and  action  taken  codes.  [Ref.  6:  pp.  17-19] 

1.  Malfunction  Description  Codes 

The  three  digit  numeric  malfunction  description  codes  describe,  as  accurately 
as  possible,  the  probable  cause  of  a  system  component  failure.  Malfunction 
description  codes  used  when  no  failure  was  detected  arc  displayed  in  Table  3. 
Maintenance  records  which  contained  these  codes  were  eliminated  from  the  final  data 
set.  Conditional  malfunction  description  codes  identify  component  failures  caused  by 
factors  other  than  natural  wear-out.  Examples  of  this  type  of  failure  include:  failed  due 
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to  weather,  foreign  object  damage  (FOD),  improper  or  faulty  maintenance,  etc.  To 
reduce  the  complexity  of  parameter  estimation,  conditional  failures  were  treated  as  if 
the  failure  had  occurred  through  normal  usage.  The  components  studied  had  a  low 
percentage  of  conditional  failures.  In  fact,  less  than  three  percent  of  all  observed 
failures  were  conditional. 


TABLE  3 

EXCLUDED  MALFUNCTION  DESCRIPTION  CODES 


CODE  DESCRIPTION 

799  No  defect 

800  No  defect  -  removed  and'or  reinstalled  to  facilitate  maintenance 

801  No  defect  -  removed  for  modification 

803  No  defect  -  removed  for  time  change 

804  No  defect  -  removed  for  scheduled ‘maintenance 

805  No  defect  -  removed  for  pool  stock 

806  No  defect  -  removed  as  part  of  a  matched  svstem 

807  No  defect  r  removal  directed  bv  higher  authority 

811  No  defect  *  removed  during  troubleshooting 


2.  Transaction  Codes 

Transaction  codes  are  two  digit  numeric  codes  used  to  identify  the  type  of 
data  reported  by  the  maintenance  activity.  Transaction  codes  used  in  the  maintenance 

c 

reporting  system  are  listed  in  Appendix  A.  Table  4  iists  those  transaction  codes 
indicative  of  maintenance  actions  that  required  removal,  installation;  and  or  repair  of  a 
defective  component.  Maintenance  records  associated  with  these  codes  were  included  in 
the  final  data  set.  All  other  maintenance  records  were  deleted. 

3.  Action  Taken  Codes 

Action  taken  codes  are  one  digit  alphanumeric  codes  which  describe  the 
maintenance  action  performed  on  the  system,  component  or  item.  A  compiete  list  of 
action  taken  codes  is  shown  in  Appendix  B.  Action  taken  codes  associated  with  the 
removal,  installation,  and  'or  repair  of  a  defective  component  were  included  in  the  final 
data  set  and  are  shown  in  Table  5. 
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TABLE  4 

TRANSACTION  CODES  USED  IN  STUDY 

CODE 

DESCRIPTION 

12 

On-equipment  work,  including  engines,  involving  non- 
repairable  components/items  documented  as  failed  parts. 

21 

Will  be  used  when  a  repairable  component  is  removed, 
excluding  engine  components,  for  processing  at  an  I  or  D 
level  maintenance  activity  (This  code  is  used  when 
only  the  removal  must  be  documented  and  a  replacement 
is  hot  required). 

23 

Removal  and  replacement  of  a  defective  or  suspected 
defective  component  from  an  end  item,  excluding  engines 
at  the  O-level.  Additionally  this  code  is  used  for  the“ 
removal  and  replacement  of  a  complete  engine  assembly 
lor  a  defect,  suspected  defect,  or  scneduled'maintenance 
requirement. 

24 

Will  be  used  when  a  repairable  engine  component 
is  removed  for  processing  at  an  I  or  D  levei  activity 
(This  code  is  used  when  cnlv  the  removal  must  be ' 
documented  and  a  replacement  is  not  required). 

25 

_ 

Removal  and  replacement  of  a  defective  or  suspected 
defective  repairable  component  from  an  engine. 

E.  FINAL  DATA  SET 
1.  Problems 

Several  problems  were  encountered  during  the  data  reduction  process  that 
require  discussion.  First,  it  was  not  always  clear  which  flight  evolution  caused  a 
component  to  fail.  The  cause  for  this  confusion  is  the  poor  record  keeping  associated 
with  the  actual  time  of  component  failures.  For  example,  the  maintenance  data 
reporting  system  does  not  record  the  actual  time  of  component  failure.  Instead,  the 
data  is  arranged  according  to  the  time  the  maintenance  activity  received  the  component 
for  repair.  In  most,  but  not  all  cases,  the  component  was  received  by  the  maintenance 
activity  immediately  following  the  flight  causing  the  component  failure. 

Second,  several  records  indicated  multiple  failures  of  a  single  component.  For 
example,  a  single  component  may  experience  several  unrelated  malfunctions  during  a 
single  flight  evolution.  To  expedite  the  repair  and  or  replacement  of  the  foiled 
component,  the  maintenance  reporting  system  requires  the  documentation  of  each 
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TABLE  5 

ACTION  TAKEN  CODES 

CODE 

DESCRIPTION 

1 

BCM  -  maintenance  activity  not  authorized  to  conduct 
repair 

2 

BCM  -  Lack  of  equipment,  tools,  or  facilities 

3 

BCM  -  Lack  of  technical  skills  necessarv  to  complete 
repair 

4 

BCM  -  Lack  of  parts 

5 

BCM  -  Fails  check  and  test,  and  maintenance  is  allowed 
to  conduct  check  and  test  only 

6 

BCM  -  Lack  of  technical  data 

7 

BCM  -  Beyond  the  authorized  repair  depth 

8 

BCM  -  Administrative 

9 

BCM  -  Condemned,  repair  not  feasible 

C 

This  code  is  used  when  a  repairable  item  of  material 
identified  by  a  work  unit  code  is  repaired. 

F 

Failure  of  components,  items  undergoing  test. 

P 

This  code  is  entered  when  an  item  of  material  is  . 
removed  and  only  the  removal  is  to  be  accounted  for. 

Q 

This  code  is  entered  when  an  item  is  installed  and 
only  the  installation  action  is  to  be  accounted  for. 

R 

This  code  is  entered  when  an  item  of  material  is 
removed  due  to  suspected  malfunction  and  the  same  or  a 
like  item  is  reinstalled. 

*  BCM:  Action  taken  codes  labelled  1-9  are  used  for  repairable 
items  which  have  been  administratively  and  or  technically  screened 
and  tound  to  be  beyond  the  capability  of  the  maintenance  activity. 

malfunction  (i.e.,  multiple  failures  of  a  single  component).  Previous  studies  have 
treated  this  problem  as  independent  single  component  failures  with  identical 
interoccurrence  times.  The  alternative  is  to  treat  the  multiple  failure  as  a  single 
incident. 

Third,  the  underlying  assumption  of  the  exponential  distribution  is  that 
components  fail  as  a  function  of  (light  hours  flown.  This  would  imply  that  a 
component  is  operated  continuously,  which  may  or  may  not  be  the  case. 
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Fourth,  it  was  difficult  to  determine  if  a  component  was  replaced  immediately 
following  the  flight  evolution  causing  the  failure.  Additionally,  there  was  some  concern 
pertaining  to  periods  of  flight  operations  during  which  a  component  might  be  missing. 

Fifth,  the  condition  of  the  component  (i.e.,  new’  or  used)  at  the  start  of 
deployment  was  not  always  know  n.  A  related  problem  involves  the  use  of  the  same 
failure  rate  distribution  for  both  repaired  and  new  components. 

2.  Simplifying  Assumptions 

Before  proceeding  with  any  data  analysis,  it  is  necessary  to  make  assumptions 
which  can  simplify  or  eliminate  the  problems  discussed  in  the  previous  section.  These 
assumptions  are  presented  in  the  same  order  as  the  problem  areas  they  are  intended  to 
simplify. 

First,  all  failures  are  associated  with  the  sortie  during  which  they  occurred  for  ■ 
When  Discovered  Codes  A,  B,  C,  and  D.  For  ail  other  codes,  it  is  assumed  thar  the 
training  and  experience  of  the  flight  and  ground  crews  ensure  the  discovery  of  the 
failure  at  the  earliest  possible  opportunity  (i.e.,  the  failure  occurred  during  the  flight 
prior  to  discovery  ).  A  complete  listing  of  When  Discovered  Codes  is  found  in 
Appendix  C. 

Second,  all  multiple  failures  of  a  single  component  are  treated  as  a  single 
component  failure.  It  is  believed  that  multiple  failures  would  affect  the  turnaround 
time  required  to  repair  the  component  and  may  cause  delays  within  the  repair  pipeline. 
However,  this  should  not  affect  the  rate  at  which  the  component  fails. 

Third,  it  is  assumed  that  cornponents  operate  continuously,  at  least  in  a 
standby  mode.  If  this  is  not  the  case,  the  sortie-dependent  models  should  provide  a 
better  estimate  of  expected  failures. 

Fourth,  the  study  assumes  that  components  are  replaced  immediately  and  that 
gaps  in  coverage  are  not  allowed.  In  fact,  it  is  extremeiy  unlikely  that  an  aircraft 
would  operate  without  any  of  the  components  listed  in  this  thesis. 

Fifth,  all  components  are  considered  new  at  the  start  of  the  deployment.  For 
the  Exponential  and  Geometric  models  this  assumption  is  not  necessary  due  to  the 
"memoryiess  property"  they  exhibit.  Additionally,  the  impact  on  the  Wcibull 
distribution  should  be  minimal  if  the  number  of  observed  failures  is  large. 

Finally,  to  reduce  the  complexity  of  the  problem,  it  is  necessary  to  assume 
that  repaired  components  exhibit  the  same  failure  characteristics  as  new  components 
(i.e.,  have  the  same  failure  rate). 
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3.  Format 


a.  Parameter  Estimation 

To  obtain  the  best  possible  estimate  of  the  time  between  failures,  the 
observations  within  the  initial  data  base  were  sorted  by  aircraft  bureau  number,  then 
Julian  date  and  the  event  time.  Failure  times  were  measured  between  consecutive 
component  failures  on  a  specific  aircraft.  These  measurements  were  never  made 
between  component  failures  on  two  different  aircraft.  The  information  included  in  the 
final  data  set  is  the  aggregate  totals  and  includes  a  censoring  indicator,  number  of 
sorties  flown,  flight  time  without  failure,  and  flight  time  with  failure. 

A  censor  indicator  equal  to  zero  indicates  that  a  component  failure  was 
observed.  The  data  entries  associated  with  this  type  of  observation  include  the  number 
of  flights  flown  inclusive  of  the  event  causing  the  failure,  total  flight  time  flown  since 
last  failure  (does  not  include  flight  time  of  event  causing  failure),  and  the  flight  time  of 
the  event  causing  the  failure. 

A  censor  indicator  equal  to  one  represents  an  event  which  is  right  hand 
censored.  Right  hand  censoring  refers  to  an  observation  of  a  component  which  does 
not  fail  prior  to  the  termination  of  the  deployment  (i.e.,  the  data  set  documents  a 
component  s  survival  time,  but  provides  no  information  pertaining  to  the  actual  time  of 
component  failure).  The  data  entries  associated  with  this  type  of  observation  arc  the 
total  number  of  flights  flown  and  the  total  flight  time  since  the*  last  component  failure. 

b.  Model  Comparison 

The  intent  of  this  thesis  is  to  determine  each  model's  ability  to  predict 
failures  for  the  components  listed,  and  compare  those  predictions  with  the  actual 
component  failures  observed  during  the  deployment  of  the  squadrons  aboard  the 
LSS  KENNEDY  and  USS  INDEPENDENCE.  To  accurately  assess  these  capabilities, 
it  was  necessary  to  identify  component  failures  as  they  occurred  during  the  deployment. 
Thus,  the  data  set  used  to  provide  the  graphical  and  quantitative  comparison  of  the 
predicted  vs  actual  component  failures  was  different  than  the  data  set  used  to  estimate 
the  model  parameters.  The  latter  data  set  was  sorted  by  Julian  date  and  event  time 
only.  The  accumulated  number  of  sorties  flown  and  flight  houi  totals  were  recorded 
with  the  observation  of  each  failure,  such  that,  comparisons  could  be  made  with  the 
expected  failures  generated  by  the  failure  rate  models. 


4.  Summary  Statistics 

Tabic  6  summarizes  the  two  data  sets  and  provides  such  information  as  the 
number  of  aircraft  assigned,  sorties  flown,  flight  hour  totals,  component  failures,  etc. 


TABLE  6 

SUMMARY  STATISTICS 


DATA  SET 

AIRCRAFT  SORTIES 

FLT  HRS  SORTIE 
FLT  HRS  MEAN  STD. DEV. 

MODEL 

30  3489 

7608.S  2.18 

.98 

VALIDATION 

44  3012 

7030.4  2.33 

r 

.54 

wuc 

MODEL  SET 
FAILURES 

.  VALIDATION  SET 
FAILURES 

734H100 

105 

58 

46X1600 

34 

49 

4622100 

49 

.  6 

56X2500 

85 

79 

632Z100 

49 

57 

74A1500 

107 

75 

74A5MOO 

118 

104 

5772200 

34 

33 

6918100 

81 

93 

MODEL  SET  =  VF-l  14  and  VF-213 
VALIDATION  SET  =  VF-11,  VF-14,  VF-31,  and  VF-32 
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III.  FLIGHT  HOUR  MODELS 


A.  OVERVIEW 

The  two  models  discussed  in  this  chapter  describe  component  failures  as  a 
function  of  flight  hours  only.  The  first  flight  hour  model  to  be  examined  is  the 
exponential,  which  is  the  model  currently  used  by  the  Navy  to  describe  the  distribution 
of  the  time  between  component  failures.  This  model  has  some  weaknesses  inherited 
from  the  properties  of  the  exponential  distribution.  The  effect  these  problem  areas 
might  have  on  the  model's  ability  to  accurately  predict  component  failures  will  be 
discussed  in  the  next  section.  The  three  models  offered  as  alternatives  should  provide 
improvement  to  some  or  all  of  these  problem  areas,  and  the  comparisons  made  in 
Chapter  V  should  shed  new  light  on  the  magnitude  of  these  deficiencies. 

The  second  flight  hour  model  discussed  is  the  Weibull.  It  is  offered  as  an 
alternative  method  of  describing  the  time  between  failures  when  flight  hours  are  the 
only  factors  contributing  to  the  component  malfunction. 


B.  EXPONENTIAL  MODEL 
1.  Weaknesses 

There  are  two  major  weaknesses  associated  with  this  model  that  cause  some 
concern.  First,  the  exponential  is  a  flight  hour  model  and  as  such  it  is  insensitive  to  the 
effects  of  the  sortie  on  the  component  s  time  of  failure.  In  Chapter  V  it  will  be  shown 
that-  the  expected  number  of  component  failures  for  the  time  interval  (0,t]  is 
proportional  to  t,  so  that  the  expected  number  of  failures  for  a  specific  flight  hour  total 
is  not  affected  by  the  average  sortie  length  (i.e.,  the  number  of  sorties  flown). 

Second,  the  exponential  distributions  are  "NO  WEAR"  distributions,  which 
indicates  that  a  component's  probability  of  surviving  a  flight  of  duration  t  is  the  same 
for  new  components  as  it  is  for  components  which  are  "S"  hours  old  (i.e.,  the 
conditional  survival  probability  is  equal  to  the  unconditional  probability  of  survival). 
This  is  the  "memory  less"  property  of  the  exponential  distribution,  and  implies  that  the 
prior  use  of  a  component  has  no  effect  on  the  time  of  its  failure. 

The  same  concern  is  reflected  in  the  exponential  distribution's  constant  hazard 
rate  r(t)  (i.e.,  rfu-X  for  all  values  of  t).  As  such,  the  failure  rate  X  would  not  be 
aflected  by  the  time  on  tost  t.  It  is  believed  that  certain  components  exhibit  hazard 
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rates  that  are  decreasing  functions  of  time  (t).  In  fact,  many  pieces  of  electronic 
equipment  are  routinely  exposed  to  bum-in  periods  to  make  the  component  less 
susceptible  to  failure  when  put  into  operation.  Conversely,  it  can  be  argued  that  other 
components  have  hazard  rates  that  increase  with  time  (i.e.,  a  component  is  more  likely 
to  fail  as  it  experiences  some  wear  and  tear).  The  exponential  distribution  is  not 
capable  of  modelling  these  differences. 

2.  Maximum  Likelihood  Results 

The  maximum  likelihood  estimate,  X',  of  the  exponential  failure  rate 
parameter  X  (failures/hour)  is  found  by  iteratively  solving  equation  3.1:2 

^  ^  texpf-Xt;)  Tj.tj^O,  X>0,  m<n 

>  T.-'  J - — 1 —  ‘I  (eqn  3.1) 

j  1  Y  l*exP(-Xt)  i=  1,2,. ...n  j-  1,2 . m 

where,  the  left  hand  side  of  equation  3.1  equals  the  aggregate  flight  time  for  those 
sorties  that  did  not  experience  a  specific  component  failure,  tj  is  the  duration  of  the 
sortie  causing  the  failure,  n  is  the  total  number  of  records,  censored  and  uncensored, 
within  the  final  data  set,  and  rn  is  the  number  of  uncensored  observations  (failures,. 

If  the  sortie  lengths,  tj,  are  equal  to  t  for  all  j,  then  equation  3.1  can  be 
simplified  to  the  closed  form  expression  given  in  equation  3.2.  1  hus,  the  average  sortie 
length  can  be  used  to  provide  a  "rough"  estimate  of  the  maximum  likelihood  estimator 
X: 

,  -ln(FTWF.'TFT) 

X  — -  (eqn  3.2) 

t 

where,  FTWF  (Flight  Time  Without  Failure)  equals  the  sum  of  the  Tj's  for  i=  l,...n, 
TFT  equals  the  total  flight  time  for  the  deployment,  and  t  equals  the  average  sortie 
length. 

The  maximum  likelihood  estimate  for  X'  was  found  using  the  FORTRAN 
program  listed  in  Appendix  E.  The  parameter  estimates  and  95  percent  confidence 
limits  for  these  estimates  are  displayed  in  Table  7. 


‘For  the  derivation  of  this  maximum  likelihood  equation,  the  reader  is  referred  to 
Appendix  D. 


TABLE  7 

EXPONENTIAL  MAXIMUM  LIKELIHOOD  RESULTS 

WUC 

X' 

ffx' 

734H100 

.0140 

.0018 

(.0106,-0175) 

46X1600 

,0048 

.0010 

(.0029, .0067) 

4622100 

.0065 

.0011 

(.0043..0086) 

56X2500 

.0114 

.0021 

(.0073,  .0155) 

632Z100 

.0067 

.0012 

(.0043..0090) 

74A1500 

.0143 

.0019 

(.0105, .0181) 

74A5M00 

.0159 

.0024 

(.01 13, .0206) 

5772200 

.0048 

.0011 

(.0026,-0070) 

6918100 

.0108 

.0015 

(.0079,-0137) 

***  CL=* 

95  PERCENT  CONFIDENCE  LIMITS 

3.  Quantile-Quantile  Plots 

A  graphical  method  often  used  to  determine  an  empirical  data's  goodness  of 
fit  with  respect  to  a  theoretical  distribution  is  the  quantile-quantile  plot.  The  technique 
requires  the  bbservations,  X^Xj.-.-.X  ,  to  be  ordered  from  smallest  to  largest. 
Empirical  quantiles  are  then  computed  for  these  order  statisiics  based  on  a  formula, 
such  as  Qi  =  (i-.5j  n,  where  Qj  is  the  empirical  quantile  for  the  X^  order  statistic  and  n 
equals  the  total  number  of  observations.  These  empirical  quantiles  are  plotted  against 
the  theoretical  quantiles  to  produce  a  graphical  estimate  of  goodness  of  fit. 

Some  adjustments  to  the  above  technique  must  be  made  to  overcome  the 
problems  of  censoring.  The  first  adjustment  relates  to  the  problem  associated  with 
incomplete  knowledge  of  the  actual  time  oi  component  failure.  In  Chapter  II,  it  was 
mentioned  that  the  failure  occurred  during  the  interval  T;  to  T  +  t.  To  simplify  this 
problem,  it  is  assumed  that  the  midpoint  of  that  interval  represents  the  actual  time  of 
failure. 


The  second  adjustment  that  needs  to  be  made  accounts  for  those  observations 
that  are  right  hand  censored.  The  Kaplan* Meier  method  for  quantile  estimation  will 
account  for  those  observations  and  will  provide  a  more  accurate  estimate  for  the 
empirical  quantiles.  This  method  also  requires  that  each  observation,  censored  and 
unceni.  ;:ed.  be  ordered  from  smallest  to  largest.  The  empirical  quantiles  (Q;)  are 
computed  using  equation  3.3: 

Qj*  l-[(n+.5)/n]  II  [(n-k  +  .5)/(n-k+ 1.5)]  for  i  in  I  (eqr.  3.3) 

k  in  J 
k£i 

where,  I  is  the  set  of  all  unccnsored  observations,  and  h  is  the  total  number  of 
observations,  censored  and  uncensored.  [Ref.  7:  p.  234] 

An  exponential  quantile-quantile  plot  was  constructed  for  each  of  the  nine 
components  studied.  The  quantile-quantile  plots  for  work  unit  codes  734H100  and 
56X2500  are  displayed  in  Figures  3.1  and  3.2,  and  the  remaining  quantile-quantile  plots 
are  presented  in  Appendix  F.  The  y-axis  of  the  plot  represents  the  quantiles,  F(T^  t), 
and  the  x-axis  represents  the  time  between  component  failures.  The  asterisks  identify 
the  uncensored  observations,  the  *  s  identify  the  locations  of  the  censored 
observations,  and  the  theoretical  distribution  evaluated  using  k'  is  defined  by  the  solid 
line.  The  empirical  data  is  fit  well  by  the  theoretical  distribution  if  the  data  points  are 
tightly  packed  around  the  solid  line.  Data  points  which  deviate  from  the  line  can 
provide  information  about  areas  of  weakness.  • 

Figure  3.1  illustrates  the  best  exponential  fit  of  all  the  quantile-quantile  plots, 
but  the  plot  still  shows  the  lower  to  mid  quantiles  associated  with  the  empirical  data 
points  overestimate  the  quantiles  computed  using  the  theoretical  distribution.  Figure 
3.2  provides  a  better  representation  of  the  other  quantile-quantile  plots.  In  even-  plot, 
the  empirical  quantiles  overestimate  the  theoretical  quantiles  for  quantiles  less  than  .5, 
and  underestimate  the  theoretical  quantile  for  quantiles  greater  than  .75. 

Figure  3.3  displays  the  Weibull  cumulative  distribution  function  for  a  equal  to 
.5,  1,  and  2.  The  reader  should  be  aware  that  a  Weibull  distribution  with  shape 
parameter  (a)  equai  to  1  is  the  exponential  distribution.  Figure  3.3  shows  that  for 
values  of  time  less  than  t0,  the  cumulative  distribution  for  the  curve  represented  by 
a  -.5  is  greater  than  the  exponential  curve  («=  l),  and  for  values  of  time  greater  than 
tQ  it  is  less  than  the  cumulative  distribution  of  the  exponential  curve.  This  is  the 
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identical  situation  displayed  in  each  of  the  nine  quantile-quantile  plots.  This  would 
suggest  that  the  distribution  representing  the  time  between  component  failures  might 
be  better  modelled  with  a  Weibull  distribution  with  parameters,  X>  0  and  0 5a  <  1. 


Figute  3.3  Weibull  Cumulative  Distribution  Function. 


C.  WEIBULL  MODEL 
1.  Properties 

The  next  flight  hour  model  that  will  be  discussed  is  the  Weibull  distribution, 
which  can  be  a  viable  alternative  to  the  exponential  distribution  since  it  is  capable  of 
modelling  many  different  failure  rates.  The  major  difference  between  the  Weibull  and 
Exponential  models  is  that  the  former  is  usually  not  "memorylcss which  indicates  that 
the  component  s  previous  usage  can  be  exploited  in  the  determination  of  the  expected 
failure  time. 
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The  two  parameters  associated  with  the  Weibull  distribution  are  the  scale 
parameter  (X)  and  the  shape  parameter  (a),  with  each  shape  parameter  defining  a 
specific  failure  rate  distribution.  As  such,  the  Weibull  family  of  distributions  has  the 
capability  to  model  many  different  types  of  component  wear.  For  example,  a  values 
greater  than  one  describe  a  family  of  "WEAR  OUT"  distributions  (increasing  failure 
rate  functions);  an  a  value  equal  to  one  describe  the  "NO  WEAR"  exponential;  and  a 
values  less  than  one  describe  a  family  of  "WEAR  IN"  distributions  (decreasing  failure 
rate  functions).  Consequently,  the  Weibull  distribution  provides  an  additional 
modelling  flexibility  that  should  reduce  the  latter  two  areas  of  concern  associated  with 
the  exponential  distribution. 

The  W'eibull  cumulative  distribution  function  (Figure  3.3)  and  the  exponential 
quantile-quantile  plots  indicate  that  a  Weibull  model  with  scale  parameter  (X>0)  and 
shape  parameter  (a<l)  might  provide  appropriate  descriptors  of  the  failure  rate 
distribution  for  each  of  the  components  studied.  As  mentioned  above,  a  shape 
parameter  less  than  one  is  indicative  of  a  "WEAR  IN"  family  of  distributions. 
Distributions  of  this  type  are  characterized  by  a  hazard  rate  function,  r(t),  which  is 
decreasing  in  t,  (i.e.,  the  component  s  failure  rate  decreases  as  t  increases).  Since,  the 
majority  of  components  studied  are  electronic  components  (AVIONICS),  this 
assumption  of  a  decreasing  failure  rate  is  probably  an  improvement  over  the  constant 
failure  rate  of  the  exponential  model. 

2.  Maximum  Likelihood  Results 

Maximum  likelihood  solutions  can  be  obtained  by  setting  the  partials  equal  to 
zero,  and  solving  for  X'  and  a'  using  nonlinear  optimization  techniques  similar  to  the 
one  described  in  Appendix  G.  The  partial  derivatives  of  the  log-likelihood  function 
with  respect  to  the  parameters  X  and  a  are  displayed  beiow:3 


dK  aXq~1((Ti  +  ti)q-Tic]ST(t1) 

^x”2-  i-s^tj) 

i 


-XaX^Tj0 

j 


dK 

da 


{(Xdj  +  ti)jalog[X(Ti  4-  ti)l-(XT-)qlog(  XTjJ  ST  (tj) 

dst^V 


-  XO-Tj)alog(XTj) 
i 


3The  derivation  of  these  equations  is  given  in  Appendix  0. 
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where,  a,X>0;  T^tj.TjiO;  i“l,2,...,m  (the  set  of  uncensored  observations  -  failures); 
and  j“m+  l,m+2,...,n  (the  set  of  censored  observations). 

The  maximum  likelihood  estimates  for  X'  and  a'  were  found  using  the 
FORTRAN  program  listed  in  Appendix  G.  The  parameter  estimates  and  95  percent 
confidence  limits  for  these  estimates  are  displayed  in  Table  8. 


TABLE  8 

WEIBULL  MAXIMUM  LIKELIHOOD  RESULTS 

wuc 

X' 

CLX 

a 

CLa 

734H100 

.0149 

(.0084, .02 16) 

.8079 

(.8014, .8146) 

46X1600 

.0044 

(.0006, .0082) 

.6664 

(.6627,.6703) 

4622100 

.0066 

(.0024, .0108) 

.7708 

(.7666, .7750) 

56X2500 

.0119 

(.0063,  .01 76) 

.8170 

(.8114..8227) 

632Z100 

.0067 

(.0026, .0109) 

.7970 

(.7928, .SOU) 

74A1500 

•  .0165 

(.0084, .0246) 

.7066 

(.6985, .7147) 

74A5M00 

.0164 

(.0104, .0224)  ' 

.8844 

(.8784,.8904) 

■  5772200 

.0045 

(.0009, .0082) 

.7385 

(.7349,.7422) 

6918100 

.0112 

(.0058..0167) 

.8142 

(.8087..S196) 

***  CL  = 

95  PERCENT  CONFIDENCE  LIMITS 

3.  Quantile-Quantile  Plots 

The  parameter  estimates  listed  in  Table  8  were  used  to  construct  Kaplan- 
Meier  quantile-quantile  plots  for  each  of  the  nine  components  studied.  Two  of  these 
quantile-quantile  plots  are  displayed  in  Figures  3.4  and  3.5,  and  the  other  seven  are 
included  in  Appendix  H.  Both  Figures  graphically  illustrate  that  the  data  is  fit  well  by 
the  Weibull  distribution,  with  the  worst  case  lit  represented  by  Figure  3.4.  This  plot 
shows  that  the  Weibull  distribution  has  a  tendency  to  underestimate  the  empirical 
quantiles  below  ten  percent,  but  fits  the  data  extremely  well  for  quantiles  greater  than 
ten  percent. 


The  results  of  the  Weibull  quantile-quantile  plots  are  not  conclusive,  but  offer 
credibility  to  the  use  of  the  Weibull  model  as  an  alternative  method  of  modelling 
component  failure  rates  when  flight  time  is  the  only  contributing  factor. 
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Figure  3.5  Weibuii  Quantile-Quantile  Plot:  WLC  56X2500. 
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IV.  SORTIE  DEPENDENT  MODELS 


A.  OVERVIEW 

The  two  models  that  will  be  discussed  in  this  chapter  differ  from  the  models  in 
Chapter  III  in  that  component  failures  can  be  affected  by  stresses  incurred  as  a  result 
of  the  sortie.4  These  stresses  refer  to  events  that  occur  during  each  flight,  regardless  of 
the  sortie  length,  and  would  include  events  such  as  equipment  initialization,  takeoffs, 
landings,  etc.  The  first  sortie-dependent  model  to  be  discussed  is  the  Geometric  model, 
which  determines  a  constant  probability  that  a  component  fails  on  a  given  sortie. 

The  second  sortie-dependent  model  that  will  be  briefly  looked  at  is  the  Mixed 
model  described  in  reference  5.  This  model  could  also  be  classified  as  a  flight,  hour 
model  since  it  determines  the  probability  of  a  component  failure  based  on  both  sortie 
and  flight  hour  contributing  factors.  It  is  incorporated  in  this  section  because  it  also 
utilizes  the  geometric  distribution  to  determine  the  probability  that  the  first  failure 
occurs  on  the  sortie.  The  major  difference  between  these  two  models  is  that  the 
Mixed  model's  probability  of  failure  for  a  given  sortie  is  not  constant.  In  fact,  it  is  a 
function  of  a  constant  failure  probability  caused  by  sortic-rclatcd  stress,  and  a 
probability  of  failure  which  is  dependent  on  the  sortie  length. 


B.  GEOMETRIC  MODEL 
1.  Properties 

The  basic  premise  of  the  Geometric  model  is  that  the  component  failure 
results  due  to  sortie-related  stress  and  is  insensitive  to  differences  in  sortie  length.  This 
assumption  can  then  be  used  to  formulate  the  componer.'s  survival  or  failure  for  a 
given  sortie  as  a  Bernoulli  trial,  where,  the  component  fails  with  probability  p  or 
survives  with  probability  1-p.  The  geometric  distribution  is  used  to  determine  the 
probability  of  observing  the  first  failure  on  the  s1*1  Bernoulli  trial. 

The  geometric  distribution  is  in  many  ways  similar  to  the  exponential 
distribution  and  might  be  considered  the  exponential's  discrete  analog.  Like  its 
continuous  counterpart,  the  geometric  distribution  is  rncmorvless"  and  has  a  constant 
failure  rate,  which  would  indicate  that  the  Geometric  model  also  ignores  the  history  of 


VA  sortie  refers  to  one  complete  flight  evolution,  take-off  to  full-stop  landing. 


component  usage.  Any  improvement  ottered  by  this  model  would  be  attributed  to  the 
effects  of  the  sortie  on  the  component's  failure. 

It  should  be  apparent  that  a  low  variability  in  sortie  length  (i.e.,  the  sortie 
length  is  approximately  constant)  vvould  imply  that  the  exponential  distribution  is 
equivalent  to  a  geometric  distribution  with  p*  1-expi-Xt);  where  t  equals  the  average 
sortie  length.  In  this  case,  both  models  would  provide  similar  estimates  of  component 
failures.  The  flight  data  used  to  generate  the  parameter  estimates  was  acquired  from 
aircraft  flying  during  Flex-deck,  operations,  and  should  have  greater  variability  than 
would  be  expected  from  aircraft  flying  under  normal  cyclic  operations.5 

2.  Maximum  Likelihood  Results 

The  maximum  likelihood  expression  for  p',  the  estimate  of  the  constant  sortie 
failure  probability  is  given  in  equation  4. 1:6 


s-=  1,2,..  i  =  1,2, ...,n 


(eqn  4.1) 


where,  m  represents  the  number  of  uncensored  observations  (sorties  with  component 

failures),  and  the  summation  of  the  .Sj’s  for  i  ~  1,2 . n  represents  the  total  number  of 

sorties  flown.  Thus,  the  maximum  likelihood  estimator  p'  equals  the  number  oi  flights 
experiencing  a  component  failure  divided  by  the  total  number  of  sorties  flown.  This  is 
the  same,  maximum  likelihood  estimate  which  would  have  resulted  without  censoring. 
The  maximum  likelihood  estimates  p'  for  p  and  r.ne  associated  95  percent  confidence 
intervals  are  listed  in  Table  9. 

C.  MIXED  MODEL 
1.  Properties 

This  section  will  provide  a  brief  description  of  the  Mixed  model  proposed  by 
Levy  (Ref.  5].  This  particular  model  is  structured  such  that  the  sortie  and  the  sortie 
length  (flight  hours)  can  combine  to  cause  a  component  failure.  It  is  obvious  that 
factors  contributing  to  one  component's  failure  may  or  may  not  be  detrimental  to 
another  component's  operating  life.  Most  likely,  "the  life  of  electronic  components  is 
probably  affected  more  by  the  number  of  power  surges  to  which  they  are  subjected 


^During  l  lcx-deck  operations,  aircraft  sortie  ieniuh  is  not  constrained  to  a  one  to 
two  hour  launch  and  recovery  cycle. 

6  l  he  derivation  of  this  equation  is  given  in  Appendix  D. 
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than  operating  hours;  landing  gear  life  should  be  influenced  primarily  by  takeoffs  and 
landings;  and  engine  life  should  be  affected  by  rotational  speed  changes,  which  are 
accompanied  by  heat  transients  and  pressure  changes"  [Ref.  8:  p.  2J.  The  Mixed  model 
appears  to  have  an  increased  flexibilty  to  model  some  of  these  sortie  and  flight  hour 
factors. 

As  in  the  previous  section,  a  constant  probability  of  failure  is  employed  to 
model  the  effects  of  sortie  induced  stress.  In  addition,  the  Mixed  model  incorporates  a 
failure  probability  for  those  components  whose  lifetime  is  affected  by  the  flight  hours 
flown.  Combining  the  constant  sortie  failure  probability  with  the  flight  hour  probability 
results  in  the  expression  for  the  probability  of  failure  pj  for  the  i1*1  sortie,  which  is 
expressed  as  the  complement  of  a  component  surviving  the  stresses  of  the  sortie  and 
the  flight  hours  flown: 

Pi  “  l-exp(-Xtj)(l*p0)  0£po£  1,  X>0,  tj2:0  (eqn4.2) 

where,  the  flight  hour  failure  distribution  is  exponential  and  pQ  equals  the  constant 
sortie  failure  probability. 

2.  Maximum  Likelihood  Results 

The  partial  derivative  of  the  log-likelihood  function  with  respect  to  the 
parameters  X  and  pQ  are  listed  below:7 

dK_  (l-p0)tiexp(-Xti)  t 
e\  f  i-exp(-Xti)(  i-p0)  *  1  |  i 

X"  expf-Xt))  s,;-l  v L 

dpQ  Y  l-cxpt-XtiK  l-p0)  1-P0  “  1-P0 

where,  0£pQil;  X>0;  Tj,t|,Tj^0;  Sj»l,2,...;  i»  1,2 . m  (the  set  of  uncensored 

observations  *  failures);  and  j  *  m -f-  l,m  +  2,...,n  (the  set  of  censored  observations). 

Maximum  likelihood  solutions  for  X'  and  p0'  were  obtained  using  the 
FORTRAN  program  listed  in  Appendix  G.  For  each  of  the  nine  components  studied, 
the  log-likelihood  function  was  maximized  when  X'  equaled  zero,  and  pQ-  equaled  p', 
the  maximum  likelihood  estimate  for  the  constant  failure  rate  p  of  the  Geometric 


(eqn  4.3) 

(eqn  4.4) 


7The  derivation  cr  these  equations  is  given  in  Appendix  D. 
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model.  The  results  for  X'  are  included  in  Table  9.  The  reader  should  note  that  the 
values  for  p'  and  p0'  are  identical,  and  the  failure  for  X  is  zero  for  every  component. 
Since  both  models  provide  the  same  maximum  likelihood  results,  they  will  be  jointly 
referred  to  as  the  Geometric  model  when  compared  with  the  other  two  models. 


TABLE  9 

MAXIMUM  LIKELIHOOD  RESULTS 
SORTIE  DEPENDENT  MODELS 

wuc 

p'** 

CLP 

v  o 

X'* 

CLX 

734H100 

.0301 

(.0188,  .0414) 

.0000 

(.0000, .0055) 

46X1600 

.0097 

(.0038, .0156) 

.0000 

(.0000,-0027) 

4622100 

.0140 

(.0069,0211) 

.0000 

(.0000, .0034) 

56X2500 

.0244 

(.0147, .0341) 

.0000 

(.0000..0047) 

632Z100 

.0140 

(.0069, .02 11) 

.0000 

(.0000, .0033) 

74 A 1 500 

.0307 

(.0191,0423) 

.0000 

(.0000, .0056) 

74A5MOO 

.0338 

(.0228, .0448)  * 

.0000 

(.0000, .0053) 

5772200 

.0097 

(.0037, .0157) 

.0000 

(.0000, .0028;  ■ 

6918100  .0232  (.0137, .0327) 

CL-  95  PERCENT  CONFIDENCE  LIMITS 

*  MIXED  MODEL  ONLY 

**  p0'-p‘ 

.0000 

(.0000, .0045) 

V.  MODEL  COMPARISONS 


A.  PURPOSE 

The  purpose  of  this  chapter  is  to  compare  the  Exponential  model  with  the 
Weibull  and  Geometric  models  to  determine  if  the  distribution  currently  used  provides 
the  most  accurate  description  of  the  failure  process.  As  a  method  of  comparison,  the 
Chi-square  goodness-of-fit  test  was  not  considered  for  two  reasons.  First,  the  failure 
rate  information  provided  by  the  censored  observations  would  have  been  lost  because 
cf  the  impossibility  of  determining  the  exact  time  of  component  failure.  Second,  it  was 
considered  more  important  to  examine  each  model's  capability  to’ predict  the  number 
of  component  failures  resulting  from  a  specified  number  of  sorties  and  flight  hours 
flown.  The  next  section  describes  the  methods  of  estimating  failures  for  each  model. 
These  estimates  are  then  compared  with  the  actual  number  of  failures  to  determine 
each  model's  goodness-of-fit. 

B.  COMPONENT  FAILURE  ESTIMATION 

1.  Geometric  Model 

The  determination  of  the  expected  number  of  component  failures  associated 
with  the  Geometric  model  is  straightforward.  The  parameter  p  estimated  using  the 
Geometric  model  defines  the  probability  that  a  component  would  fail  during  a  given 
sortie.  Furthermore,  if  it  is  assumed  that  the  outcome  of  a  sortie  is  independent  of  the 
outcome  of  the  other  N-l  sorties,  and  p  is  constant  for  each  sortie,  then  the  probability 
law  defining  the  number  of  component  failures  for  the  deployment  is  the  sum  of  N 
independent  Bernoulli  trials.  This  random  sum  has  a  Binomial  distribution  with  an 
expected  value  equal  to  N  *  p.  Thus,  using  the  Geometric  model  the  expected  number 
of  component  failures  is  calculated  as  the  product  of  the  number  of  sorties  flown  (N) 
and  the  probability  that  the  component  fails  during  the  sortie. 

2.  Exponential  Model 

To  compute  the  expected  number  of  component  failures  using  the  Exponential 
model,  each  failure  can  be  thought  of  as  an  occurrence  within  a  Renewal  process;  a 
stochastic  process  which  counts  the  number  of  occurrences  within  an  interval  (0,t|.  It  is 
assumed  that,  the  times  between  successive  occurrences  are  independently  and 
identically  distributed  non-negative  random  variables.  As  such,  it  can  be  shown  that 
the  expected  number  of  occurrences  within  an  interval  (0,t]  is  equal  to: 
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M(t)-£Fk(t)  k- 1,2,3,... 
k 

where,  M(t)  equals  the  expected  number  of  occurrences  (failures)  within  the  interval 
(0,tj,  and  Fk(t)  is  the  k-fold  convolution  of  the  cumulative  distribution  function.  The 
derivative  of  M(t)  with  respect  to  t  defines  the  renewal  rate  of  the  process  as: 


M'(t)-m(t)«£  fk(t)  k- 1,2,3,... 
k 


Since  the  k-fold  convolution  of  the  exponential  probability  density  function  is  known 
to  be  Gamma,  m(t)  is  equal  to  the  infinite  sum  of  Gamma  probability  density  functions 
with  scale  parameter  X  and  shape  parameter  k. 


m(t)  =  £ 
k 


kkxk-le-Xx 

T(k) 


(Xx)^1 

Too" 


k-  1,2,3,... 


In  the  equation  above,  Y  ((Xx^'V'I’fk)]  for  k  =  1,2,3,...  is  easily  recognized  as  the 
Taylor's  series  expansion  of  eAx.  Thus,  the  renewal  rate  for  the  counting  process 
defined  by  exponential  interoccurrence  times  is  X.  The  expected  number  of  failures 
M(t)  is  found  by  integrating  m(t)  from  0  to  t. 

t 

\l(t)=  f  X  ds-Xt 
0 


Therefore,  the  expected  number  of  component  failures  is  equal  to  the  product  of  the 
failure  rate  X  and  the  time  on  test  t.  This  particular  Renewal  process  is  referred  to  as  a 
Poisson  process. 

3.  Weibuli  Model 

Estimating  the  number  of  expected  failures  associated  with  the  Weibuli  model 
is  not  as  easy  as  the  other  methods  discussed  above.  The  observation  of  component 
failures  is  still  a  Renewal  process,  but  the  interoccurrence  times  are  distributed 
according  to  a  more  complex  distribution  function.  Finding  the  k-fold  convolution  of 
either  the  cumulative  distribution  function  or  the  probability  density  function  is  an 
extremely  difficult  if  not  impossible  task.  For  processes  defined  by  continuous 
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interevent  times  with  finite  means  and  variances,  the  Elementary  Renewal  Theorem  a 

states  that  M(t)  is  approximately  Normal  with  mean®  t/p  and  variance®  t<T^/p^  as  t  1 

approaches  infinity  [Ref.  9:  p.  290].  Unfortunately,  this  theorem  does  not  help  estimate  3 

failures  when  t  is  small.  As  an  alternate  method  to  construct  M(t)  when  t  is  small,  a 
simple  computer  simulation  is  used. 

The  inputs  required  for  the  simulation  included  the  number  of  aircraft 
assigned  to  the  carrier  air  wings,  the  Weibull  shape  and  scale  parameters,  and  an  N 

approximate  flight  hour  distribution.  This  information  was  used  to  repeatedly  simulate  * 

the  air  wing  operations  during  the  deployment.  | 

The  aircraft  were  chosen  at  random  to  fly  a  fictitious  mission  of  some  sortie  1 

length  t.  For  the  modelling  set,  there  were  30  aircraft  assigned.  The  model  assumed  II 

that  each  aircraft  had  an  equally  likely  (1/30  percent)  probability  of  being  chosen  to  fly  . 
a  mission.  Realistically,  aircraft  are  not  chosen  at  random  to  fly  missions,  and  most  ■ 

likely  those  aircraft  recently  flown  are  better  prepared  for  future  flights.  A  weighting 

. 

scheme  could  have  been  used  to  increase  the  probabilities  of  those  aircraft  flown 

recently,  and  decrease  the  probability  fox  those  aircraft  that  were  undergoing  simulated  ( 

maintenance  activities.  3 

* 

The  duration  of  the  mission  was  determined  by  a  random  variable  generated  j 

from  a  Normal  distribution  with  the  mean  and  variance  estimated  by  the  Model  data  j 

set  (see  Table  6).  Figure  5.1  depicts  a  histogram  of  the  sortie  lengths  for  the  Model  set  3 

with  an  overlay  of  the  appropriate  Normal  distribution.  Between  the  values  of  0  to  6  r 

hours,  the  Normal  distribution  fits  the  histogram  adequately  enough  to  justify  its  use  jj 

as  the  sortie  length’s  parent  distribution.  Values  of  randomly-generated  sortie  lengths  ■ 

less  than  .5  hours  or  greater  than  6  hours  were  truncated  to  .5  hours  and  6  hours, 
respectively.  i| 

After  the  selection  of  aircraft,  the  conditional  probability  of  failure  was  ) 

computed  for  a  flight  of  sortie  length  t.  The  only  information  required  for  this  jj 

computation  was  the  component's  Weibull  shape  and  scale  parameters,  and  the  \ 

aircraft's  flight  history  since  the  last  component,  failure.  The  conditional  probability  of 
failure  could  then  be  calculated  as  described  in  Appendix  D.  The  value  for  the  2 

conditional  probability  of  failure  was  then  compared  against  a  Uniform  (0,1)  random  | 

variable,  and  if  the  uniformly-generated  outcome  was  less  than  or  equal  to  the  1 

conditional  probability  of  failure,  a  component  failure  was  simulated.  J 


45 
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Figure  5. 1  Sortie  Length-  Histogram  vs  Normal  Distribution. 

The  number  of  component  failures  was  recorded  at  500  hour  intervals  until 
the  completion  of  the  deployment.  The  deployment  was  simulated  100  times,  and  the 
average  number  of  component  failures  was  computed  for  each  500  hour  interval.  These 
averages  were  the  approximations  to  the  expected  number  of  component  failures  \1(t) 
ai  500,1000, 1500,...  hours  etc. 

C.  COMPARISONS 
1.  Methodology 

The  predicted  failures  were  compared  with  the  actual  failures  as  they  occurred 
in  each  of  the  two  deployments  iL'SS  KENNEDY  and  USS  INDEPENDENCE).  At 
this  point,  it  is  important  to  clarify  the  specific  calculations  of  the  expected  failures.  As 
described  earlier,  the  accumulated  flight  time  and  total  sorties  flown  were  recorded  with 
the  observation  of  each  failure.  Therefore,  the  exponential  and  geometric  expected 
failures  could  he  computed  using  the  formulas  listed  :n  this  chapter.  Unfortunately,  the 
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expected  number  of  Weibull  component  failures  was  not  always  defined  at  ihe 
accumulated  flight  time  of  the  actual  failure.  To  provide  such  an  estimate,  it  was 
necessary  to  interpolate8  between  the  Weibull  expected  failures  associated  with  the 
simulated  failure  times  above  and  below  the  actual  one  of  component  failure.  The 
difference  between  the  expected  and  actual  failures  was  computed  at  the  time  of  each 
component  failure.  These  differences  were  then  summed  and  divided  by  the  total 
number  of  component  failures  to  provide  an  estimate  of  the  average  prediction  error. 
An  average  difference  close  to  zero,  with  small  standard  deviation,  indicates  that  the 
model  is  a  good  predictor  of  actual  failures.  A  negative  average  shows  that  the  model 
has  a  tendency  to  underestimate  the  number  of  component  failures,  and  a  positive 
average  is  indicative  of  overestimation.  Tables  10  and  11  provide  the  average 
differences  and  corresponding  standard  deviations  for  the  "model"  and  "validation"  data 
sets  respectively. 

The  model  providing  the  best  fit  for  a  specific  component  was  chosen  based 
on  the  average  difference  between  the  predicted  failures  and  the  actual  failures.  Table 
12  provides  a  summary  listing  of  the  optimal  distributions  for  each  component  in  both 
data  sets.  The  reader  is  cautioned  that  the  optimal  was  chosen  based  on  an  average 
difference  closest  to  zero,  and  in  some  cases  the  expected  number  of  failures  was 
adequately  described  by  other  distributions  as  well  as  the  one  listed  as  optimal.  If  two 
or  more  models  provided  similar  estimates,  the  One  with  the  smallest  standard  deviation 
was  chosen  as  optimal.  The  actual  failures  were  plotted  against  the  predicted  failures  to 
further  illustrate  a  distribution's  ability  or  disability  to  estimate  failures.  These  graphs 
are  included  in  Appendix  I. 

2.  Results 

Before  proceeding  with  the  comparisons,  a  caveat  should  be  placed  on  the 
results.  The  ability  of  the  models  to  predict  the  failures  observed  during  the  deployment 
of  the  "validation"  data  set  is  contingent  on  the  accuracy  of  the  model  parameter 
estimates.  For  a  valid  comparison,  it  is  necessary  that  both  deployments  come  from  the 
same  population.  That  is  to  say,  the  flight  operations  were  conducted  under  similar 
conditions.  The  reader  should  be  aware  that  the  "model"  set  aircraft  and  the 
"validation"  set  aircraft  operated  during  different  seasons  and  in  different  geographic 
locations.  The  parameters  estimated  from  the  "model"  set  may  or  may  not  provide  the 

A  quadratic  interpolation  of  the  form  aFT4'  +  bFT-t-c=  EF  was  used.  Where.  FT 
equaled  accumulated  flight  time  ol  failure,  EF  equaled  expected  failures,  and  a.b.c  were 
parameters  estimated  using  the  known  Weibull  times  and  expected  failures. 
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TABLE  10 

MODEL  SET  AVERAGE  DIFFERENCES 


WUC«734H100 
Average : 

Exponential 

Geometric 

Weibull 

-0.31 

1.51 

3.20 

St. Dev. : 

2.58 

3.53 

2.89 

WUC (46X1600 

Average i 

0.94 

0.50 

3.65 

St. Dev. i 

1.91 

1.27 

1.93 

WUC (4622100 

Average : 

-2.87 

-1.91 

0.56 

St. Dev. s 

2.54 

1.97 

2.52 

WUC (56X2500 

Average : 

-1.37 

0.17 

1.43 

St. Dev. : 

1.94 

1.18 

1.76 

WUC(6322100 

Average  t 

-0.04 

0.08 

2.74 

St. Dev. a 

1.85 

2.12 

2.29 

WUC (74A1500 
Average  s 

-1.75 

-0.03 

5.78 

St. Dev. : 

4.77 

3.92 

5.69 

WUC (74A5M00 

Average « 

-0.53 

1.05 

0.09 

St. Dev. ( 

3.81 

4.45 

4.14 

WUC (5772200 
Average ( 

-1.49 

-1.84 

1.95 

St.Dev. : 

2.13 

1.29 

2.44 

WUC (6918100 

Average ( 

-4.36 

-2.97 

-1.75 

St.Dev. ( 

3.23 

2.53 

3.03 

true  estimate  of  the  population  parameters  required  to  accurately  predict  failures  within 
the  "validation"  data  set. 

Within  the  "model"  set,  all  of  the  models  predict  failures  adequately.  This  is 
not  unexpected  because  the  parameters  defining  the  distribution  were  estimated  using 
the  flight  data  from  the  "model"  set.  Table  1 1  shows  that  this  is  not  true  for  the 
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TABLE  11 

VALIDATION  SET  AVERAGE  DIFFERENCES 


WUC:734H1Q0 
Average i 

Exponential 

Geometric 

Weibull 

20.76 

16.44 

26.22 

St. Dev. i 

11.01 

8.85 

12.01 

WUC .46X1600 
Average  t 

-5.31 

-8.04 

-0.37 

St. Dev. t 

4.04 

5.32 

4.17 

WUC.4622100 

Average  t 

16.20 

14.57 

20.90 

St. Dev. s 

15.33 

13.95 

16.86 

WUC:  56X2500 

Average : 

-4.23 

-7.41 

0.48 

Sc. Dev. s 

4.01 

5.15 

3.49 

WUC  :632Z100 

Average  : 

-5.88 

-8.43 

-1.65 

St. Dev. j 

3.14 

4.45 

2.42 

WUC :74A1500 
Average  : 

5.82 

2.00 

17.12 

St. Dev. s 

6.10 

4.26 

8.06 

WUC  :74A5MOO 

Average  s 

8.67 

2.83 

TO. 39 

.St.Dev.  : 

3.41 

3.75 

3.74 

WUC: 5772200 

Average : 

-4.61 

-6.33 

-0.07 

St. Dev. : 

1.93 

2.40 

2.02 

WUC -.6918100 

Average : 

-2.01 

-5.81 

2.57 

St. Dev. « 

7.12 

8.43 

7.27 

deployment  defined  by  the  "validation"  data  set.  In  most  cases,  a  single  optimal  model 
can  be  identified. 

Table  10  illustrates  an  extremely  minor  difference  between  the  exponential  and 
geometric  estimates.  In  fact,  it  could  be  argued  the  geometric  and  the  exponential  wore 
interchangeable  as  predictors  of  failures  for  the  "model"  set  deployment.  For  the 
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TABLE  12 

OPTIMAL  DISTRIBUTIONS 


wuc 

MODEL  SET 

VALIDATION  SET 

734H100 

Exponential 

Geometric 

46X1600 

Geometric 

Weibull 

4622100 

Weibull 

Geometric 

56X2500 

Geometric 

Weibull 

632Z100 

Exponential 

Weibull 

74A1500 

Geometric 

Geometric 

74A5M00 

Geometric 

Geometric 

5772200 

Geometric 

Weibull 

6918100 

Weibull 

Exponential 

Inertial  Measuring-  Unit  (WUC734H100),  the  optimal  distribution  is  listed  as 
exponential  for  the  "model"  set  and  geometric  for  the  "validation"  set,  but  it  could  have 
been  geometric  in  both  cases.  During  the  deployment  represented  by  the  "validation" 
data  set  the  average  sortie  length  increased  by  approximately  8.5  percent.  As 
mentioned  earlier,  the  exponential  distribution  is  insensitive  to  changes  in  the  average 
sortie  length  and  would  provide  an  estimate  based  on  flight  hours  only.  On  th£  other 
hand,  if  the  geometric  distribution  was  correct  the  longer  sortie  lengths  would  imply 
that  a  specific  flight  hour  total  was  accomplished  with  fewer  sorties  flown.  It  would 
provide  a  smaller  estimate  than  would  the  Exponential  model.  For  Work  Unit  Code 
734H100,  the  Exponential  model  was  unable  to  account  for  the  longer  average  sortie 
length  (in  the  "validation"  set)  and  overestimated  the  actual  number  of  failures. 

The  Geometric  model  appears  to  provide  the  best  fit  for  components  that  are 
not  necessarily  operated  continuously.  Examples  of  components  of  this  type  wouid 
include  the  AWG-9  Radar  Transmitter  (WUC:74A  1500)  and  Detail  Data  Display 
(WUG74A5M00).  A  possible  exception  to  this  ruie  would  be  the  Inertia!  Measuring 
Unit  (WUC734H100).  Despite  its  continuous  use,  the  Inertial  Measuring  Unit's  gyros 
and  accelerometers  are  extremely  sensitive  to  the  movements  and  impacts  associated 
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with  carrier  takeoffs  and  landings,  and  it  is  extremely  likely  that  component  failures  are 
caused  by  sortie-related  stress  vice  continuous  flight  operations. 

For  continuous  use  components  within  the  “validation"  data  set,  the  Weibull 
distribution  provided  the  best  estimate  of  component  failures.  The  exponential 
distribution  was  the  optimal  predictor  for  only  one  component  (WUC:6918100).  and  in 
that  case,  the  Weibull  distribution  was  nearly  as  accurate.  From  a  readiness 
standpoint,  the  Weibull  model  might  be  considered  a  better  choice  because  of  its 
tendency  to  overestimate  the  actual  number  of  failures.  Amazingly,  the  exponential 
distribution  predicted  failures  best  for  less  than  10  percent  of  the  components  within 
the  “validation"  data  set.  In  support  of  the  current  model,  it  should  be  noted  that  the 
exponential's  estimate  of  failures  was  more  than  adequate  for  the  "model"  set. 

As  mentioned  previously,  the  Weibull  distribution  with  shape  parameter  equal 
to  one  is  exponential.  Figure  5.2  demonstrates  the  capability  of  the  Weibull  simulation 
to  provide  estimates  sinular  to  those  generated  for  the  Exponential  model.  This  figure 
show's  the  Weibull  and  exponential  predictions  for  the  Detail  Data  Display 
(WUC:74A5MOO)  with  estimated  Weibull  shape  parameter  equal  to  .91. 

It  should  not  be  surprising  that  in  some  cases,  none  of  the  models  will  be 
adequate  estimators  of  the  actual  failures.  This  type  of  phenomenon  was  observed  for 
the  Fuel  Tank  Release  Mechanism  (WUC4622100).  For  two  successive  deployments, 
total  failures  went  from  49  to  6  for  a  similar  number  of  flight  hours  and  semes  llown. 
Based  on  the  available  data,  it  is  hard  to  determine  if  49  was  too  high  or  6  was  too 
low.  This  type  of  outcome  emphasizes  the  difficulty  in  predicting  the  actual  failures  for 
an  upcoming  deployment  based  on  parameters  estimated  from  one  deployment  rather 
than  the  entire  historical  data  base. 
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Figure  5.2  Flight  Hour  Models:  Model  Set  WL'C:74A5M00. 
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VI.  SUMMARY 


A.  ADVANTAGES  AND  DISADVANTAGES 

This  section  provides  a  brief  summary  of  the  primary  advantages  and 
disadvantages  associated  with  each  of  the  four  models  discussed  in  Chapters  III  and 
IV.  The  major  advantage  of  the  Exponential  and  Geometric  models  is  that  the  number 
of  expected  failures  is  easily  calculated  using  closed  form  expressions.  Both  models 
express  failures  as  a  linear  function  of  only  one  factor,  flight  hours  (Exponential)  or 
sorties  flown  (Geometric).  There  is  a  tendency  then  to  concentrate  on  aggregate  flight 
totals  while  ignoring  the  specific  flight  history  of  each  aircraft.  As  a  result,  neither 
model  is  able  to  adjust  its  failure  estimates  to  compensate  for  differences  in  flight 
programs. 

The  Mixed  model,  on  the  other  hand,  attempts  to  account  for  failures  resulting 
from  sortie  stress  and.  or  continuous  use.  Unlike  the  previous  models,  the  expected 
number  of  component  failures  is  a  function  of  a  specific  flight  progiam.  The 
probability  of  failure  associated  with  the  Mixed  model  is  a  function  of  the  sortie 
duration.  As  such,  this  model  should  be  more  sensitive  to  differences  in  flight 
programs  and  would  provide  a  more  accurate  estimate  of  failures  -than  the  Exponential 
and  Geometric  models.  Use  of  the  Exponential  distribution  to  model  continuous-use 
components  has  weakened  the  forecasting  capabilities  of  the  Mixed  model.  Lack  of 
consideration  for  flight  hours  flown  since  the  last  component  failure  oversimplifies  the 
problem.  As  stated  previously,  the  small  variability  in  sortie  length  combined  with  the 
"memoryless"  property  of  the  Exponential  distribution  would  imply  that  the  probability 
of  flight  hour  induced  failure  is  approximately  constant  for  all  sorties  f  lown.  Since  the 
stress  related  failure  probability  is  also  constant,  this  model  would  not  necessarily 
predict  any  better  than  the  Geometric  model.  This  type  of  phenomenon  was  observed 
in  this  study. 

The  Weibull's  strongest  feature  is  its  ability  to  model  failure  rates  lor  components 
which  have  increasing,  decreasing,  or  constant  failure  rates.  Since  it  does  not  possess 
the  memoryless  property  of  the  Exponential,  Geometric  models,  the  expected  number  of 
failures  is  dependent  upon  previous  flight  histories.  Like  the  Mixed  model,  this  should 
provide  an  improved  estimate  of  component  failures.  One  disadvantage  is  a  lack  of 
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simple  mathematical  formulas  with  which  to  calculate  the  expected  failures.  However, 
in  today's  world  of  computers  this  should  not  be  a  major  problem. 

B.  CONCLUSIONS  AND  RECOMMENDATIONS 

It  is  extremely  difficult  to  estimate  the  number  of  component  failures  which 
occur  during  a  deployment  based  solely  on  the  data  provided  by  another  deployment. 
Such  an  approach  would  not  account  for  differences  in  flight  and  maintenance 
personnel,  weather  conditions,  geographical  location,  mission  type,  aircraft's 
maintenance  history,  etc.  For  example,  the  number  of  actual  failures  observed  for  the 
Fuel  Tank  Release  Mechanism  went  from  forty-nine  to  six  for  a  similar  number  of 
flight  hours  and  sorties  flown.  To  minimize  the  effects  caused  by  such  variations  in 
deployments  it  is  recommended  that  these  models  be  evaluated  using  an  enlarged  data 
base. 

Surprisingly,  in  this  study,  the  model  currently  used  by  the  Navy  (Exponential) 
did  not  provide  the  best  estimate  of  component  failures.  The  Geometric  provided  a 
better  fit  for  components  which  were  not  operated  continuously,  and  the  Weibull 
performed  better  when  the  components  were  operated  continuously.  In  fact,  overall, 
the  Exponential  was  the  least  effective  model  for  these  nine  components. 

The  results  obtained  using  the  Weibull  simulation  were  encouraging.  The  model 
produced  satisfactory  estimates  with  a  minimal  number  of  assumptions,  data  input,  and 
software  coding.  While  the  results  obtained  using  this  model  were  not  overwhelming, 
an  improvement  in  prediction  capability  was  observed. 

Two  possible  directions  for  improving  the  current  method  of  estimating  demand 
(i.c.,  component  failures)  are  the  construction  of  a  universal  model  and  the  use  of  high 
resolution  simulations.  The  first  alternative  suggests  the  use  of  a  model  which  cun 
describe  components  with  different  failure  rate  functions,  and  incorporates  a 
methodology  utilizing  sortie-induced  stress  and  flight  hours  flown  to  compute  an 
estimate  of  component  failures.  The  model  currently  used  does  not  have  this 
capability. 

A  slight  modification  to  the  Mixed  model  described  in  Chapter  IV  might  provide 
the  modelling  flexibility  desired.  A  failure  rate  distribution,  such  as  the  Weibull.  could 
be  used  to  describe  continuous-use  failures  instead  of  the  Exponential.  Such  a  change 
would  ensure  that  the  component’s  failure  probability  could  be  influenced  by  prior  use. 
That  is  to  say.  the  probability  of  failure  lor  a  given  sortie  would  now  be  a  function  of 


the  sortie  length  and  the  aggregate  flight  time  since  the  last  component  failure.  The 
strengths  of  the  Weibull  distribution  combined  with  the  constant  probability  of  failure 
(sortie)  associated  with  the  Mixed  model  should  provide  an  estimate  of  demand  which 
is  more  sensitive  to  variations  within  and  between  flight  programs. 

The  price  paid  for  this  improvement  in  forecasting  ability  is  an  increase  in  the 
complexity  of  the  expected  demand  computations.  With  the  abundance  of  computers 
within  the  fleet  this  should  not  be  a  problem.  Perhaps,  ten  to  twenty  years  ago  it  was 
necessary  to  rely  on  simple  models  for  prediction  and  the  cnide  estimates  they 
provided.  This  is  not  the  case  today;  computer  simulations  can  be  used  when  simple 
mathematical  formulas  are  not  available.  As  described  in  Chapter  V,  a  low  resolution 
simulation  was  used  to  forecast  the  component  failures  associated  with  the  Weibull 
renewal  process.  The  advantage  of  such  simulations  is  the  ability  to  include  other 
factors  which  may  affect  a  component's  failure  rate.  Examples  of  these  factors  would 
include  number  and  type  of  landings,  mission  type,  maintenance  programs,  weather 
conditions,  etc.  It  is  safe  to  say,  that  these  more  complicated,  high  resolution  models 
could  include  more  information  about  the  factors  which  cause  components  to  fail,  and 
should  ultimately  result  in  better  forecasts  of  demand. 
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APPENDIX  A 
TRANSACTION  CODES 


CODE  DESCRIPTION 

00  Used  to  report  an  inventory  gain. 

02  Used  to  report  q  change  in  the  readiness  reportable 
status  of  an  equipment. 

03  Used  to  report  an  equipment  loss. 

1 1  a.  On-equipment  work  not  involving  a  removal  or 
replacement  of  a  defective  component  or  item. 

b.  On  supproting  engine  documents  not  having  a  removal 
of  a  defective  or  suspected  defective  component  or  item 
when  the  engine  is  not  specifically  identified  to  a 
specific  aircraft. 

c.  Used  at  the  O  or  I  level  maintenance  activities  when 
closing  out  a  maintenance  action. 

12  On-equipment  work,  including  engines,  involving  non¬ 
repayable  components  items  documented  as  faired  parts. 

14  Removal  of  a  nondefective  component,  item,  excluding 
cannibalization,  from  an  engine  to  be  processed  at  the 
O-level  maintenance  activity. 

15  Installation  of  a  nondefective  component/item, 
excluding  cannibalization,  on  an  engine  to  be  processed 
at  an  O-Tevel  maintenance  activity. 

16  Removal  of  a  nondefective  component  item,  excluding 
engine  components  items  and  a  cannibalization  to  be 
processed  at  an  O-level  maintenance  activity. 

17  Installation  of  a  nondefective  component/item 
excluding  engine  component  items  and  cannibalization. 

18  a.  The  removal  and  replacement  of  nondefective 
engines  and  component,  items  to  accomplish  a 
cannibalization  action. 

b.  The  removal  and  replacement  of  those  consumeable 
components  items  subject  to  a  scheduled  removal  interval 
or  items  ol  supplv  significance. 

c.  The  removal  and  replacement  of  a  nondefective  engine 
component  for  cannibalization  at  the  I  -level  only. 

19  a.  The  removal  and  replacement  of  nondefective 

engine  component  to  accomplish  a  cannibalization  at  the 
O-Tevel  onlv. 

b.  The  removal  and  replacement  of  those  consumeable 
engine  components  subject  to  a  scheduled  removal  interval 
or  items  ol  supply  significance. 

21  Will  be  used  when  a  repairable  component  is  removed. 

excluding  engine  components,  lor  processing  at  an  I  or  D 
level  maintenance  activity  (This  code  is  used  when 
only  the  removal  must  be  documented  and  a  replacement 
is  not  required). 

33  Removal  and  replacement  of  a  defective  or  suspected 

defective  component  from  an  end  item,  excluding  engines 
at  the  O-level. 

Additionally,  this  code  is  used  for  the  removal  and 
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replacement  of  a  complete  engine  assembly  for  a  defect, 
suspected  defect,  or  scheduled  maintenance  requirement. 

24  Will  be  used  when  a  repairable  engine  component 
is  removed  for  processing  at  an  I  or  D  level  activity 
(This  code  is  used  when  only  the  removal  must  be 
documented  and  a  replacement  is  not  required). 

25  Removal  and  replacement  of  a  defective  or  suspected 
defective  repairaole  component  from  an  engine. 

30  Used  to  document  components  processed  through  the 
I -level  maintenance  activity  for  check,  test,  and 
service. 

31  Work  performed  on  a  removed  repairable  component  with 
no  failed  parts  or  awaiting  parts  documented  in  the 
Failed/ Required  Material  mocks  of  the  VIDS.MAF. 

32  Work  performed  on  a  removed  repairable  component  with 
failed  pans,  awaiting  parts,  or  cannibalization 

actions  documented  in  the  Failed/Required  Material 
blocks. 

39  Close  out  fer  man  hours  or  awaiting  parts  at  the 
I -level  maintenance  activity. 

41  Technical  directive  compliance  with  no  part  number 
change. 

47  Technical  directive  compliance  with  a  part  number 
change. 

72  Will  be  useed  to  report  SCIR  data  bv  the  reporting 

custodian  when  transient  maintenance  is  performed  bv 
other  than  the  reporting  custodian. 
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APPENDIX  B 
ACTION  TAKEN  CODES 


CODE  DESCRIPTION 

1  BCM  -  maintenance  activity  not  authorized  to  conduct 
repair 

2  BCM  -  Lack  of  equipment,  tools,  or  facilities 

3  BCM  -  Lack  of  technical  skills  necessary  to  complete 
repair 

4  BCM  -  Lack  of  parts 

5  BCM  -  Fails  check  and  test,  and  maintenance  is  allowed 
to  conduct  check  and  test  only 

6  BCM  -  Lack  of  technical  data 

7  BCM  -  Beyond  the  authorized  repair  depth 

8  BCM  -  Administrative 

9  BCM  -  Condemned,  repair  not  feasible 

A  Items  of  repairable  Materia!  or  Weapon.  Support 

Svstem  Discrepencv  Checked  No  Repair  Required.  This 
code  is  used  for  all'discrepencies  which  are  checked 
and  found  that  either  the  reported  dehciencv  cannot 
be  duplicated,  or  is  operating  within  allowable 
tolerances. 

B  Repair  or  replacement  of  attaching  units,  seals. 

gaskets,  etc.,  that  are  not  integral  parts  of  work 
unit  coded  items. 

C  This  code  is  used  when  a  repairable  item  of  material 
identified  by  a  work  unit  code  is  repaired. 

D  This  code  is  used  to  closeout  a  V1DS  MAF  when 

component  repair  is  to  be  performed  at  another  facility. 

F  Failure  of  components,  items  undergoing  test. 

J  This  code  is  used  when  ap  item  is  calibrated  and  found 

serviceable  without  need  for  adjustment. 

K  This  code  is  used  when  an  item  must  be  adjusted  to 
meet  calibration  standards. 

L  This  code  is  used  when  a  maintenance  action  must 

be  stopped  or  delaved  while  awaiting  parts  which  are 
not  available  locallv,  and  a  component  goes  into 
an  awaiting  parts  status. 

N  This  code  is  used  bv  an  organizational  activity 

when  it  becomes  ne'cessarv  to  closeout  a  maintenance 
action  during  o:  at  the  en'd  of  a  reporting  period 
for  any  reason. 

P  This  code  is  entered  when  an  item  of  material  is 

removed  and  only  the  removal  is  to  be  accounted  for. 
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Q  This  code  is  entered  wh?n  ap  item  is  installed  and 
only  the  installation  action  is  to  be  accounted  for. 

R  This  code  is  entered  when  an  item  of  material  is 

rerngved  due  to  suspected  malfunction  and  the  same  or  a 
like  item  is  reinstalled. 

S  This  code  is  entered  when  an  item  of  material  is 

removed  to  facilitate  other  maintenance  and  the  same 
item  is  reinstalled. 

T  This  code  is  used  when  an  item  of  material  is 

removed  and  replaced  for  a  cannibalization  action. 

Y  This  code  is  u$ed  when  the  time  expended  in  locating 
a  discrepency  is  great  enough  to  warrant  separating 
the  trouolesnoot  time  from  the  repair  time. 

Z  This  code  is  used  when  actuallv  treating  corroded 

items,  and  includes  cleaning,  treatment,  pruning  and 
painting. 
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APPENDIX  C 

WHEN  DISCOVERED  CODES 

CODE  DESCRIPTION 

A  Before  flight  -  mission  aborted  -  failure  discovered  by  aircrew 

B  Before  flight  -  no  abort  -  failure  discovered  by  aircrew 

C  Inflight  -  mission  aborted 

D  Inflight  -  no  abort 

E  After  flight;  between  flight  -  failure  discovered  by  aircrew 

F  Pilot/NFO  weekly  inspection 

G  Acceptance,  transfer  inspection 

H  Between  flights  -  failure  discovered  by  ground  crew  • 

J  Daily  inspection 

K  Preflight,  daily; preflight,  postflight,  or  turnaround  inspection 

L  Special  inspection 

M  Calender  odd/major;phase  inspection 

N  Calender  even  inspection 

O  Administrative 

P  Functional  checkflight 

Q  Conditional  inspection 

R  Quality  assurance  inspection 

S  Oil  analysis  inspection 

U  Modification, SDLM  Overhaul'Airline  maintenance 

V  Related  maintenance  action 

W  In-shop  repair  disassembly  for  maintenance 

X  Test  bench.engine  test  stand  operation 

Y  Receipt  or  withdrawal  from  supply 
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APPENDIX  D 

MAXIMUM  LIKELIHOOD  DERIVATIONS 


1.  MAXIMUM  LIKELIHOOD  ESTIMATORS 
a.  Properties 

The  method  used  in  this  thesis  to  estimate  the  model  parameters  is  maximum 
likelihood.  As  such,  it  is  necessary  to  review  two  important  properties  of  this  estimator. 
First,  if  a  random  sample  of  size  n  is  taken  of  a  random  variable  X  whose  probability 
mass  function,  p(x),  depends  on  an  unknown  parameter,  0,  then  the  probability  law  for 
the  maximum  likelihood  estimator  is  asymptotically  normal  (as  n  gets  large)  with: 

MEAN  =  0  and  VARIANCE  =  (nK2)*1  (eqn  D.l) 

where,  K2=  E[((d/c0)lnp(x))2]  for  a  discrete  random  variable  X,  and  n  equals  the  total 

number  of  observations.  Thus,  the  maximum  likelihood  estimator  is  asymptotically 

unbiased;  that  is,  the  expected  value  of  the  estimator  is  equal  to  the  unknown 

♦ 

parameter  0.  Furthermore,  the  variance  of  the  estimator  is  the  Cramer- Rao  lower 
bound,  and  represents  the  smallest  possible  variance  of  an  unbiased  estimator. 
[Ref.  10:  pp.  372,379] 

The  second  property  describes  the  probability  law  for  maximum  likelihood 
estimators  when  the  X;  s  for  i=  l,...,n  are  not  identically  distributed.  In  the  models 
discussed  in  this  thesis,  the  probability  mass  function  for  each  X;  was  dependent  on 
one  or  more  other  random  variables.  Nevertheless,  the  distribution  of  the  maximum 
likelihood  estimator(s)  is  still  asymptotically  normal  with: 

MEAN  -  0  and  VARIANCE  =  (X^2)'1  (eqn  D.2) 

i 

where,  —  E[((5  .  c0)lnp(x))2]  is  evaluated  using  the  specific  parameters  defining  the 
i1^  observation.  The  diilerence  in  the  variance  in  equation  D.l  and  the  variance  in 
equation  D.2  is  caused  by  probability  mass  functions  for  the  X-  s  which  are  not 
identically  distributed.  For  a  complete  derivation  of  this  result,  the  reader  is  referred  to 
W.  J.  Heintzelman  s  1975  paper  on  Determining  the  Failure  Rate  When  Failure  Times 
Are  Sot  Known  Exactly  [Ref.  1 1]. 


b.  Confidence  Intervals 

Confidence  intervals  for  the  unknown  parameters  can  then  be  calculated  using 
the  maximum  likelihood's  probability  law  described  in  the  previous  section.  These 
limits  are  defined  as: 

CONFIDENCE  LIMITS  «  *  <T0  ) 

where,  0'  is  the  maximum  likelihood  estimator  for  0,  Oq'  is  the  standard  deviation  of 
the  maximum  likelihood  estimator,  and  Zi.a/2 ’s  l'a/2  quantile  of  the  standard 
normal  distribution. 

c.  Formulation  of  Probability  Mass  Functions 

The  generalized  format  used  to  construct  each  model's  probability  mass  function  for  X- 
is: 

p(X|)  =  [P(failure  occurred)] *'xi  x  [P( right  hand  censored)Jxi 

The  random  variable,  X;,  is  an  indicator  variable  equal  to  zero  if  the  data  entry 
contains  a  failure,  and  equal  to  one  if  the.  observation  is  right  hand  censored.9  The 
following  variables  are  used  to  derive  the  maximum  likelihood  estimate!  s)  in  each  of 
the  models: 

•  n  :  Total  number  of  records  in  the  final  data  set 

•  m  :  Total  number  of  records  containing  component  failures 

•  n-m  :  Total  number  of  records  which  were  right  hand  censored 


2.  EXPONENTIAL  MODEL 
a.  Properties 

The  probability  density  function  fit),  survival  function  S(t),  and  hazard  rate 
function  r(t)  for  the  Exponential  distribution  are  defined  below. 

fit)  =  Xe’Xt  X>0 

S(t)=Pr(TSt)=e'Xt 


X 


9 Right  hand  censored  observations  identify  those  components  that  did  not  fail 
prior  to  tne  termination  ol'the  deployment. 
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b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  outlines  the  formulation  of  the  maximum  likelihood  equation 
used  to  find  X',  the  maximum  likelihood  estimate  for  the  failure  rate  parameter  X.  The 
final  data  set  contains  information  supporting  two  record  types.  The  first  record  type 
identifies  observations  of  components  failing  during  some  specified  interval  of  time,  T- 
to  Tj  +  tj.  The  probability  equation  representative  of  this  type  of  record  is  given  in 
equation  D.3. 

P(Tj£  t^T-t-ti)  =  exp(-XTi)  *  (l-exp(-Xtj))  (eqn  D.3) 

The  second  record  type  pertains  to  an  observation  of  a  component  surviving 
past  some  time  T-.  These  censored  observations  contain  no  information  regarding  the 
actual  time  of  component  failure.  Equation  D.4  defines  the  probability  relationship 
associated  with  this  type  of  observation. 

S(T)  =  expt-XT)  (eqn  D.4) 

The  probability  mass  function,  p(Xj).  for  the  Exponential  modei  is  obtained  by 
combining  equations  D.3  and  D.4.  The  random  variable  X;  equals  zero  when  a  failure 
is  observed  and  one  when  the  data  is  right  hand  censored  (i.e.,  component  survived 
until  the  end  of  the  deployment). 

p(Xj)=  [exp(-XT)]xi  x  [exp(-XTj)  x  ( l-expl-Xr)]- 1_V  (eqn  D.5) 

where. 

•  Xj={0,1},  X>0,  and  Tj(  tiO 

•  Tj  equals  the  aggregate  fligiit  time  without  component  failure  for  the  ith  record 

•  tj  equals  the  flight  duration  of  the  sortie  during  which  the  failure  was  observed 

The  likelihood  function  for  ptv),  L(p(\j))  can  be  expressed  as  the  product  of 
the  n  probability  mass  functions.  The  maximum  likelihood  estimate  lor  X  is  the  value 
of  X'  which  maximizes  L(p(x->).  Since  the  logarithmic  function  is  a  monotonically 
increasing  function,  the  X'  which  maximizes  L(p(\j))  is  identical  to  the  X'  which 
maximizes  the  log-likelihood  function  K(p(.\j)).  The  optimal  value  for  X'  is  found  by 
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taking  the  derivative  of  the  log-likelihood  function  with  respect  to  X,  setting  it  equal  to 
zero,  and  solving  for  X. 


Likelihood  function: 

L(P(V)“P  [exp(-XTj)]xi  x  [expi-XTj)  *  (l-expl-Xtj^'V  (eqn  D.6) 
i-  1.2,3 . n 

Log-likelihood  function: 

K(p(Xi))«X  -  ^Tixi*  ^Ti(!-Xi)  +  (l-Xj)log[ I-expK-Xt;)]  {eqn  D.7) 


The  partial  derivative  of  the  log-likelihood  function  with  respect  to  X  is: 


dK  _ 


[(l-xi)tiexp(-Xti)] 

(l-expl-Xr)] 


(eqn  D.8) 


Equation.  D.9  is  the  result  of  setting  Equation  D.8  equal  to  zero  and 
rearranging  the  summations  to  eliminate  the  x^s):  The  left  hand  side  of  Equation  D.9 
is  summed  over  the  set  (i)  of  all  records,  censored  and  uncensored.  But,  the  right  hand 
side  is  only  summed  for  the  uncensored  observations  (j). 


V  x  =V  tiexP{'Xti) 
i  ‘  "j*  l-exp(-Xt.) 


Tj.tj  £0,  X>0,  m<n 
i®  1.2 . n  j  =  1,2,. ...m 


(eqn  D-9) 


c.  Maximum  Likelihood  Variance 

As  stated  in  Section  D.  1,  the  maximum  likelihood  estimator  for  X  is 
approximately  normally  distributed  with  mean  equal  to  X  and  variance  equal  to  1  £K.* 
for  i=  l,...,n.  In  this  model: 

,  (t  exp(-Xt  )  -  T.)^  ■> 

K;2  = - *— ■ - '  - x  [expl-XTj)  x  exp(  -Xtj)]  +  Ti'expi-XT;) 

1  l-exp(-Xt-)  1  11  1 

where,  is  evaluated  using  the  maximum  likelihood  estimate  of  X  (i.e.,  X  =  X  ).  The 
standard  deviation  of  the  maximum  likelihood  estimator  and  associated  confidence 
limits  displayed  in  Table  7  were  computed  using  the  equation  listed  above. 


3.  WEIBULL  MODEL 

a.  Properties 

The  probability  density  function  f^t),  survival  function  S(t),  and  hazard  rate 
function  r(t)  for  the  Weibull  distribution  are  displayed  below. 

^t)-Xaata*1exp(-(>Lt)a]  X,  a>0  t£0 

S(t)  =  exp[-(Xt)a] 

rfO-X^t®"1 

b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  describes  the  formulation  of  the  equation  used  to  find  the 
maximum  likelihood  estimates  of  the  Weibull  scale  (X)'and  shape  (a)  parameters.  The 
two  different  record  types  supporting  this  equation  are  identical  to  those  identified  in 
the  preceding  section.  The  key  difference  between  the  two  models  is  that  the  Weibull 
distribution  possesses  the  "memoryless"  property  for  a  equal  to  one  only.  As  such,  the 
conditional  survival  probability  is  usually  not  equal  to  the  unconditional  survival 
probability.  The  conditional  survival  probability  for  the  Weibull  distnbution  is  defined 
as:10 

„  ...  cxp(.|X(Ti  +  tj)|<‘) 

Ti  1  exp[-(XT)a) 

The  conditional  probability  of  failure  during  the  interval  T  to  T.  + 1.  can  be 
expressed  as  l-Sj(t-),  the  complement  of  the  conditional  survival  probability.  The 
probability  mass  function  for  the  Weibull  model  is  then  formed  by  combining  the 
conditional  probability  of  a  component  failure  within  the  interval  T-  to  T|+t-  (Xj  =  0) 
with  the  probability  that  a  component  survives  time  T-  (Xj=  1). 

p<xi)=[l-ST(ti)J1  xi  x  expl-fXT/l'i  (eqn  D.ll) 

where. 

•  Xj=  {0,1},  X>0,  and  T.  i  £0 

•  T^  equals  the  aggregate  flight  time  without  component  failure  for  the  ith  record 


iUThe  subscripted  T,  is  used  to  indicate  the  conditional  survival  probability 


j  i.  to  uauu  iiiuiv,uic  uil  uuiiuiuuiuil  oi 

associated  with  a  component  that  has  survived  Tj  hours  of  operation 
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•  tj  equals  the  flight  duration  of  the  sortie  during  which  the  failure  was  observed 


The  log-likelihood  function,  K(p(Xj)),  is: 


KfpiXj))  -  X  ( l-x.){log[  1-ST  (tj)]}  4-  ^[-(XT)0] 
i  1 


i-  1,2,3,. ...n 


(eqn  D.12) 


The  partial  derivatives  of  the  log-likelihocd  function,  KfpfXj)),  with  respect  to  the 
parameters  X  and  a  are  displayed  below.  They  have  been  simplified  by  removing  the 
xi(s)  and  summing  over  the  appropriate  records.  The  set  i=  l,2,3,...,n  represents  all  data 
records  (censored  and  uncensored),  and  the  set  j  =  m+  l,m+2,...,n  includes  only  those 
records  which  are  censored. 

— -V  o^q~  ‘t(Tj  +  tj)q-Tiq|ST  ( tj)  , 

«■  ~  I^> T) 

^  V  ![VTi-Hi)|”lo8[X(T^ti)H>.Til«»loglXTi)ST(ti) 

da  “  l-ST(t|)  "  ’  y>.Tj)  logiXTj) 

c.  Maximum  Likelihood  Variance 

The  following  equations  were  used  to  compute  the  variance  for  the  maximum 
likelihood  estimators,  X'  and  a'.  These  variances  were  then  used  to  construct  the 
confidence  limits  displayed  in  Table  8. 

''y 

with  respect  to  X': 

K,2(  V)  -  <«*«->)-’  -  (  +  T^expHXT/l) 

*’5T.'v 

with  respect  to  a': 

r  2(Tr  )_<  t^i^iogixp-,  +  ti)i-(XT1>ui0g(XTi,>isT(tl))2 


where,  p  equals  ((XT^logt ;X Ti)J-expi-(>-Ti )ttJ. 


i  nf  V\ si'Oivi'C'.Y s.  v’T'v.vw'Wt/Ct.V.  v  ’  fytataAA.*.  vW'.vWaW 


4.  GEOMETRIC  MODEL 

a.  Properties 

The  probability  mass  function,  Pr(F=s),  and  survival  function,  Pr(Fis),  for  the 
Geometric  distribution  are  given  below. 

Pr(F-  s)=  pqs"^  0  £  p  i  1,  q=  1-p  and  s=  1,2,3,... 

Pr(F^s)=qs‘1 

where,  F  is  the  random  variable  representing  the  observ  ation  of  the  first  failure,  and  q 
is  the  probability  of  the  component  surviving  the  sortie.  Equations  D.13  and  D.14 
illustrate  the  Geometric  distribution's  "memoryless"  and  constant  failure  rate  properties 
respectively. 

Pr(F>a  +  b  F>a)  =  (qa+b):qa  =  qb  =  Pr(F>b)  (eqnD.13) 

Pr(F=  s).Pr(F  £  s)  =  (pq5'1)  q^*  =  p  (eqn  D.14) 

b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  formulates  the  equation  used  to  find  the  maximum  likelihood 
stimate  of  the  constant  failure  probability  (p)  associated  with  any  given  sortie.  For 
the  Geometric  model,  the  two  event  types  which  can  occur  include  components  failing 
during  the  s,b  sortie,  or  components  failing  during  some  unknown  sortie  alter  the  stb 
sortie.  The  probability  expression  associated  with  the  first  type  of  event  is  given  in 
equation  D.  15,  and  equation  D16  defines  the  relationship  lor  those  observations  where 
the  component  survived  at  least  s  sorties. 

Pr( F  —  s>  =  pq5'1  (eqn  D.15) 

Pr(F>s»  =  qs  (eqn  D. 16) 

The  piobabilitv  mass  function  for  this  model  is  obtained  by  combining 
equations  D.15  and  D.16.  Again,  X  is  allowed  to  take  the  values  zero  and  one  only, 
‘  with  zero  indicating  a  component  failure  on  the  stb  sortie  and  a  one  indicating  a 

I  component  failure  alter  the  stb  sortie. 


.  . 

ft  fcV  afKvf.  M t  AjiAi  Aj»  +1  m  V,**P„*V 


ft  r  m 


^  t  Jt  f  A_ji  w  ^  ^  A.*  1 


p{Xj)  -  [p<  1  -p)Sf 1 J 1  *Xi  X  [(  1  -p)si]xi 


(eqn  D.17) 


where, 

•  OSpi  1,  q=  1-p,  Sj=  1,2,3,... 

•  S:  equals  the  aggregate  number  of  sorties  flown  inclusive  of  the  sortie  causing 
the  failure  when  X:  =  0  and  equals  the  aggregate  number  of  sorties  without  a 
component  failure  when  *=  1 

The  maximum  likelihood  estimate  p'  of  the  unknown  sortie  failure  rate  will  be 
found  by  taking  the  derivative  of  the  log-likelihood  function  (K(p(Xj))]  with  respect  to 
p,  setting  the  derivative  equal  to  zero  and  solving  for  p. 


K(p(Xj))-]T  (l-Xj)ln(p)  f  (Sj-l)O-Xj)ln(l-p)  +  sftlnO-p)  (eqn  D.18) 
i 

i=  1,2, 3,.  ..,n 


^K_T,(l-xi)  [(s^Hd-X;)]  ( s ;  x  | 

dp  p  (1-p)  (1-p) 

The  partial  derivative  listed  above  can  be  simplified  by  removing  the  indicator 
variables  and  adjusting  the  summations  to  produce  the  closed  form  expression  lor  p' 
given  in  equation  D.  19. 


V  Jl_ 
T  >*P 


m 

p'  *  — *=; -  i*  1,2 n  (eqn  D.19) 

i 

where,  j=  1,2,3,.. .,m  is  the  set  of  uncensored  records  (i.e. ,  component  failures),  and 
i=  1,2,3... .,n  is  the  set  of  all  records,  censored  and  uneensored.  Thus,  the  maximum 
likelihood  estimator  p'  equals  the  number  of  flights  experiencing  a  component  failure 
divided  by  the  total  number  of  sorties  flown, 
c.  Maximum  Likelihood  Variance 

The  maximum  likelihood  estimators  variance  was  computed  using  the 
formulas  described  in  Section  D.l.  The  equation  delining  the  (jeomeiric  model's 
expression  for  is  given  below: 
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(eqn  D.20) 


Kj2  *  l(l-pSi)2(l*p)V3]/p  +  Sj2(l-p)sr2 

where,  Kj2  is  evaluated  at  p  =  p'.  These  variances  were  used  to  compute  the  confidence 
limits  listed  in  Table  9. 

5.  MIXED  MODEL 

a.  Properties 

Combining  a  constant  sortie  failure  probability,  p0,  with  a  flight  hour  probability 
distribution  results  in  an  expression  for  the  probability  of  failure,  pj,  for  the  i1*1  sortie. 
This  is  expressed  as  the  complement  of  a  component  surviving  both  the  stresses  of  the 
sortie  and  the  flight  hours  flown: 

Pj  =  l-exp(-Xt-)(l-p0)  0<po  £  1,  X>0,  tj£0  (eqn  D. 21) 

where,  the  flight  hour  failure  distribution  is  Exponential. 

The  model  is  then  formulated  in  the  same  manner  as  the  Geometric  model 
described  in  the  preceding  section.  The  probability  that  the  first  failure  occurs  on  the 
sth  sortie  is  expressed  in  equation  D.22.  Note  that  the  p-^^  are  not  necessarily  equal 
fori=l . s-1 . 

Pr(F=  s)  *  ps  n.(l-pj)  i—  1.2 . s-1  (eqn  D.22) 


Equation  D.23  results  from  using  the  expression  for  p-  given  in  equation  D.21 
to  rewrite  equation  D.22. 

Pr(F  **  s)  -  [l-exp(-Xt){  l*p0)](( l-p0)s‘!exp(-XT)]  Tpt^O  (eqn  D.23) 

b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  describes  the  formulation  of  the  equation  used  to  find  the 
maximum  likelihood  estimates  for  the  failure  rate  parameter  (X)  and  the  constant  sortie 
probability  of  failure  (p0).  The  probability  relationship  supporting  the  observation  of 
the  first  failure  on  the  s1^1  sortie  is  given  in  equation  D.23,  and  the  expression  for  the 
observation  of  a  failure  at  some  unknown  time  after  s1^  sortie  is  listed  below. 
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Pr(F>  s)  -  I7(l-Pj)  -  exp(-XTXl-p0)s 


(eqn  D.24) 


The  probability  mass  function  p(Xj)  for  tnis  model  results  from  combining 
equations  D.23  and  D.24.  Notice  that  when  X  equals  zero  the  probability  mass 
function  for  the  Mixed  model  is  identical  to  the  probability  mass  function  expressed  for 
the  Geometric  model  in  equation  D.17,  and  when  p0  equals  zero  the  probability  mass 
function  is  equivalent  to  the  Exponential  probability  mass  function  given  in  equation 
D.5.  The  above  statement  implies  that  the  Mixed  model  has  the  capability  of 
modelling  the  Geometric  or  Exponential  cases  as  the  situation  would  require. 

p(Xj)=  {[l-exp(-Xti)(l-p0)][(l-p0)si*1exp(-XT1)]j  l*xi  *  [( l-p0)si  exp(-XTi)Jxi 

where, 

'•  0£-poSl,X>0;  sj-  1,2,3 . and  Tj(  tiO 

•  Tj  equals  the  aggregate  flight  time  without  component  failure  for  the  ith  record 

•  t:  equals  the  flight  duration  of  the  sortie  during  which  the  failure  was  observed 

•  S:  equals  the  aggregate  number  of  sorties  flown  inclusive  of  the  sortie  causing 
the  [ailure  when  Xi^O  and  equals  the  aggregate  number  of  sorties  without  a 
component  failure  frhen  Xj  =  If 

The  log-likelihood  function  K(p(Xj)),  and  the  partial  derivatives  of  the  log- 
likelihood  functions  with  respect  to  the  parameters  X  and  pQ  are  given  below.  The 
maximum  likelihood  estimates  for  X  and  pQ  were  obtained  utilizing  the  Quasi-Newton 
method  described  in  Appendix  G.  The  log-likelihood  function,  K(p(X;)),  is: 

K(pfxi))  =  X(l.xi)(ln[l-exp(-Xti)(l-p0)]  +  ((sl-l)ln(l-p0)]-XTi}+xi{{s1ln(l-p0)]-XTi} 
i 

i-  1.2.3 . n 


The  partial  derivatives  of  the  log-likelihood  function,  K(p(Xj)).  have  been  simplified  by 

removing  the  x^s)  and  summing  over  the  appropriate  records.  The  set  i  =  1,2.3 n 

represents  all  data  records  (censored  and  uncensored),  and  the  set  j  =  m+  l.m  +•  2 n 

includes  only  those  records  which  are  censored. 


(eqn  A. 25) 


y  (l-PoHiexpC-Xtj)  yT_ 

*mf  l-cxjX-XtiKI-Po)’  ‘  j  j 

£K_  ^  ex^-Xtj)  srl  y _ Sj_ 

^p0  j  l-exjX-XtjKl-Po)  1*P0  j  1-P0 


(eqn  A. 26) 


c.  Maximum  Likelihood  Variance 

The  equation  for  the  variance  of  the  probability  distribution  for  p0'  is 
identical  to  equation  D.20  when  X  equals  zero.  The  equation  for  the  Kj2  for  X'  is  given 
below. 


Ki2(X'>  - [(( l-p0)ti]  p^}2  x  p0(  l-p0)sr1  +  Tj2(  l-p0)s 


where,  K-2  is  evaluated  at  X  =  X'  =  0.  These  variances  were  used  to  compute  the 
confidence  limits  for  X  listed  in  Table  9. 


"I 


cywAivw^V-v.v.v. 


APPENDIX  E 

SINGLE  PARAMETER  OPTIMIZER 


*  COMPUTES  MAXIMUM  LIKELIHOOD  VALUE  FOR  A  SINGLE  UNKNOWN 

*  PARAMETER  USING  GOLDEN  RATIO  LINE  SEARCH. 


IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/RL/FTNF ,FTF, PARAM , PSD 
COMMON/ INT/ CENSOR, NF, COUNT 

DIMENSION  CENS0R(200) .NF(200) ,FTNF(200) ,FTF(2Q0) ,PARAM(2) 
INTEGER  CENSOR, COUNT, TFLAG 
REAL*8  LB,LC 


DATA  LB/ .0001D0/ITER/1/K/1/STEP/ .00001D0/ALPHA/ .OOOIDO/ 


*  INITIALIZE 
COUNT" 1 

PARAMd^.OOOOlDO 

TFLAG=0 


*  COMPUTATION  OF  SEARCH  INTERVAL  TAU.  TAU  IS  APPROXIMATELY 

*  .618  AND  TAUI  IS  APPROXIMATELY  .382 

TAU=(5**.5-l)/2 

TAU1=1-TAU 


*  INPUT  DATA 

5  READ  (3,50, END=99)  CENSOR ( COUNT ) ,NF( COUNT ), FTNF( COUNT ), FTF( COUNT) 

COUNT=COUNT+1 
GO  TO  5 

99  COUNT=COUNT-l 

*  BOUNDING  PHASE  -  SWANN'S  METHOD 

CALL  ML (LB, FI) 

6  ALPHA=ALPHA+(2**K)*ST£P 

CALL  ML ( ALPHA, F2) 

*  IF  NEW  VALUE  IS  LESS  THAN  PREVIOUS  VALUE  TERMINATE  BOUNDING 

IF  (F2.LT.F1)  GO  TO  10 
F1=F2 
K=K+1 
GO  TO  6 
10  UB=ALPHA 


*DEFINE  INITIAL  POINTS  &  ASSOCIATED  FUNCTIONAL  VALUES 


B1=TAU1*(UB-LB)+LB 
B2=TAU*(UB-LB)+LB 
CALL  ML(B1 , FI ) 

CALL  ML(B2,F2) 

*  COMPARES  THE  FUNCTIONAL  VALUES  OF  THE  INTERNAL  PTS .  B1  &  B2 ,  AND 

*  SELECTS  THE  VALUE  WITH  THE  LARGEST  FUNCTIONAL  VALUE  AS  THE  UPDATED 

*  PARAMETER  ESTIMATE 

150  IF  (F1.GT.F2)  THEN 
PARAM(2)=B1 
UB=B2 
B2=B1 
F2=F1 

Bi=TAUl*(UB-LB)+LB 
CALL  ML ( B 1 , F 1 ) 

1TER=ITER+1 
GO  TO  200 
END  IF 
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IF  (F2.GE.F1)  THEN 
PARAM(2)SB2 
LE-B1 
B1=B2 
F1=F2 

B2=TAU* ( UB-LB ) +LB 
CALL  ML(B2,F2) 

ITER=ITER+1 
GO  TO  200 
END  IF 

200  CALL  TERM (FI ,F2 , TFLAG) 

PARAM( 1 )=PARAM( 2) 

IF  (TFLAG. EQ.l)  GO  TO  275 
GO  TO  150 

*  WRITE  OUTPUT 
275  WRITE  (6,300) 

WRITE  (6,310)  PARAM(l) 

CALL  SD 

WRITE  (6,305)  PSD 
UC=PARAM(1)+1.96*PSD 
LC=PARAM ( 1 ) - 1 . 96*PSD 
WRITE  (6,320)  LC,UC 

*  FORMAT 

50  FORMAT  (II , IX, 14 , 5X, F6 . 1 , IX, F3 . 1 ) 

300  FORMAT  ( 1 OWUC : ' / ' OCONVERGENCE  CRITERION  ESTABLISHED 1 ) 

305  FORMAT  (' OSTANDARD  DEVIATION  IS  :  ',F9.8) 

310  FORMAT  ( ' OMLE  FOR  UNKNOWN  PARAMETER  IS:  1 ,F9.8) 

320  FORMAT  ('0.95  AS5YHPT0TIC  CONFIDENCE  LIMITS  ARE  : ' , F5 . 4 , '  , ' , F5 . 4) 
9999  STOP 
END 


SUBROUTINE  TERM  (FI ,  F2’,  TFLAG) 

*  DETERMINES  IF  TERMINATION  CONDITION  HAS  BEEN  SATISFIED 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/RL/FTNF , FTF , PARAM , PSD 
COMMON/ INT/ CENSOR , NF , COUNT 

DIMENSION  CENSOR(200) ,NF(200) ,FTNF(200) ,FTF(200) ,PARAM(2) 
INTEGER  TFLAG 

TEST=ABS ( F1-F2 ) 

IF  (TEST. LE. .00003)  TFLAG=1 

RETURN 

END 


SUBROUTINE  ML(PT,FVAL) 

*  COMPUTES  LOG-LIKELIHOOD  FUNCTIONAL  VALUE 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/RL/FTNF, FTF, PARAM, PSD 
COMMON/ INT/ CENSOR , NF , COUNT 

DIMENSION  CENSOR (200) ,NF(200) ,FTNF(200) ,FTF(200) ,PARAM(2) 
INTEGER  CENSOR, COUNT, TFLAG 
REAL*8  K1 

FVAL=0 

*  INPUT  FLIGHT  DATA 

DO  5  1=1, COUNT 

*  UNCENSORED  RECORDS 
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*  *  tna>  *  *-  #-  *  * 


IF  ( CENSOR ( I). EQ.O)  THEN 
Kl*EXP(-PT*FfF(I)) 
FVAL=FVAL+LOG (1-Kl) -PT*FTNF ( I ) 
GO  TO  5 
END  IF 

*  CENSORED  RECORDS 

FVAL= FVAL - PT*  FTNF ( I ) 

5  CONTINUE 

RETURN 

END 


SUBROUTINE  SD 

COMPUTES  ASSYMPTOTIC  STANDARD  DEVIATION  FOR  THE  MAXIMUM  LIKELIHOOD 
ESTIMATE 


IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/ RL / FTNF , FT  F , P ARAM , PS  D 
COMMON/ INT/CENSOR. NF , COUNT 

DIMENSION  CENS0R(200) ,NF(200) ,FTNF(200) ,FTF(2C0) ,PARAM(2) 
INTEGER  CENSOR, COUNT 
REAL*8  KISQ,K1 ,K2,K3 


SUMKI=0. 


INPUT  FLIGHT  DATA 
DO  5  1=1, COUNT 


UNCENSORED  RECORDS 

IF  (CENSOR(I) .EQ.O)  THEN 
Kl=EXP(-PARAh(l)*FTF(I) 
K2=EXP( -PARAM( 1 ) *FTNF( I 
K3s((FTF(I)*Kl/(l-Kl)J 


) 


END 


.  FTNF (!))*( ( 

KISQ=K3*{K2*(  1-Kl )  )  +  (frTNF { I )  *FTNF  ( 
GO  TO  6 
IF 


n^rKi/< 


1-K1))-FTNF(D) 


CENSORED  RECORDS 

KISQ=FTNF ( I ) *FTNF ( I ) *K2 
SUMKI=SUMKI+KISQ 
CONTINUE 


PSD=SQRT ( I/SUMKI ) 

RETURN 

END 


SUBROUTINE  SD1 


COMPUTES  ASSYMPTOTIC  STANDARD  DEVIATION  FOR  THE  MAXIMUM  LIKELIHOOD 
ESTIMATE  FOR  THE  SORTIE  DEPENDENT  MODELS 


IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/ RL/ FTNF , FTF , PARAM , PSD 
COMMON/ INT/CENSOR . NF , COUNT 

DIMENSION  CENSOR(200) ,NF(200) ,FTNF(200) ,FTF(200) ,PARAM(2) 
INTEGER  CENSOR, COUNT 
REAL*8  L,LL1 , LL2 


*  INITIALIZE  PARAMETER  ESTIMATES 
P= . 0232D0 
L= . OOOODO 
SUMKI1=0 .ODO 
SUMKI2=O.ODO 
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V.V./.V.'AV 


*  -INPUT  FLIGHT  DATA 

DO  5  1=1, COUNT 
C1=(1-P*NF (I ) ) 

C2=(1-P)**NF(I) 

C3=C2/(1-P) 

C4“C3/(l-P) 

C5=C4/(1 -P) 

C6=((l-P)*rTF(I) )/P-FTNF(I) 
SUMKI1=SUMKI1+(C1*C1*C5)/P+NF 
SUMKI2=SUMKI2+C6*C6*P*C3+FTNF 
5  CONTINUE 

*  COMPUTE  STANDARD  DEVIATIONS 

PSD1=SQRT ( 1/ SUMKI 1 ) 

PSD2=SQRT (1/ SUMKI 2 ) 

*  COMPUTE  CONFIDENCE  LIMITS 

LL1=P-1.96*PSD1 
UL1=P+1.96*PSD1 
LL2=L-1.96*PSD2 
UL2=L+1 . 96*PSD2 


Hi 


*NF ( I  )*C4 
*FTNF(I )*C2 


320 


330 


OUTPUT 

WRITE 

FORMAT 

C) 

WRITE 

FORMAT 

C.F5.4) 


(6,320)  LL1 , UL1 
('0.95  ASSYMPTO' 

(6,330)  LL2 ,UL2 

('0.95  ASSYMPTOTIC  CONFIDENCE  LIMITS  FOR  LAMBDA: 


95  ASSYMPTOTIC  CONFIDENCE  LIMITS  FOR  P :  ' , F5.4 , ' , ' ,F5.4 


,F5.4, 1  , 


RETURN 

END 


APPENDIX  F 

QUANTILE-QUANTILE  PLOTS 


figure  F.l  Exponential  Quantile-Quantile  Plot:  WL'C  40X1600, 


time  between  failures 


lujAiib*. rfi.WiriV«V»V » ifaV.  >'«V|T irfw  fhiV»  w‘«‘< ^  <«*  i\V»V.’  w*^  n'uVJvr*  JCrfiiVy'  if  >'  wvi  >■  AHbW#i  *Cr 
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Figure  F.S  W'cibull  Quantile-Quantile  Plot:  W'UC  40X1600. 
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TIME  BETWEEN  FAILURES 
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TIME  BETWEEN  FAILURES  HRS 
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Figure  F.12  Weibull  Quantile-Quantile  Plot:  WL'C  74A5MOO. 
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Figure  F.14  Weibull  Quantile-Quantile  Plot:  WUC  6918l<»0. 
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APPENDIX  G 

DUAL  PARAMETER  OPTIMIZER 


1.  QUASI-NEWTON  METHOD 

The  Quasi-Newton  method  was  used  to  determine  the  maximum  likelihood 
estimates  for  the  Weibull  and  Mixed  models.  This  method  was  derived  from  the 
Newton's  method  through  the  Taylor's  expansion  of  the  objective  function,  f(x),  where 
x  is  the  n  x  1  column  vector  containing  values  for  the  n  parameters  being  estimated. 
The  Newton  method  uses  the  quadratic  approximation  to  fix),  equation  G.l,  evaluated 
at  xW  to  iteratively  determine  an  optimal  solution. 

f{x)~f{x^)+  Vf(x^)^Ax  + .5Ax^72ftx^)Ax,  (eqn  G.l) 

The  next  point  in  the  sequence  is  found  by  forcing  the  gradient  of  the  quadratic 
approximation  (equation  0.2)  to  zero  and  sohing  for  x^4" 

7fIx)=7fIx(k))  +  72ax(k))Ax  =  0  (eqn  G.2) 

Allowing  Ax  to  equal  x^k+  ^-x^,  equation  G.2  can  be  solved  for  x^k+  ' 

x(k+  1)_  x(k)_y 2t^x(k)^- 1  (eqn  Q  j) 

The  two-parameter  optimizer  displayed  in  Appendix  G  did  not  use  the  Newton 
method  discussed  above  because  the  technique  can  be  unreliable  for  nonquadratic 
functions.  Instead,  a  more  robust  Quasi-Newton  method  was  used.  The  major 
difference  between  the  two  methods  is  that  the  Quasi-Newton  does  not  require  use  of 
the  inverse  Hessian  matrix  V2f(x)*^.  Instead,  it  uses  an  n*n  "A"  matrix  to  estimate 
the  inverse  Hessian  matrix.  The  A^  matrix  was  iteratively  updated  using  the  Broyden- 
Fletcher-Shanno  method,  and  should  be  a  lair  approximation  to  VN’(x~)  ,  when  .v'  is 
the  optimal  solution.  Substituting  A'k^  for  72l(x^kV*  in  equation  G.3  results  in  the 
Quasi-Newton  method  for  the  detet ruination  of  the  next  point  in  the  sequence.  This 
iterative  process  continues  until  the  objective  function  is  maximized.  [Ref  12:  pp. 
112-114] 
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2.  FORTRAN  PROGRAM 


*  COMPUTES  MAXIMUM  LIKELIHOOD  ESTIMATES  FOR  TWO  PARAMETER 

*  MODELS  USING  QUASI -NEWTON  METHOD. 

IMPLICIT  REAL*8  (A-H.0-2) 

COMMON/RL/ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 
COMMON/ INT/CONFLG . COUNT , TFLAG , CENSOR ,NF 
DIMENSION  CENS0R( 200 ) ,NF(200) ,DIR(2> ,FTNF(200> 

DIMENSION  FTF(200) ,GRAD(4) ,PARAM(4) 

DIMENSION  A(4) ,A1 (4) ,A2(4) , DELP(2) ,DELG(2) 

INTEGER  CENSOR, CONFLG, COUNT, TFLAG 

*  INITIALIZE  'A'  -  APPROXIMATION  FOR  INVERSE  HESSIAN 

DATA  A(l)/1 .ODO/A(2)/O.ODO/A(3)/0 ,0D0/A(4)/i . 0D0/ ITER/ 0/ IRES/ 21/ 

CONFLG=0 

DIR(1)=-I.0D0 

DIR(2)=-1.0D0 

*  INPUT  DATA 

rat  r  input 

WRITE  (6,*)  'COUNT' , COUNT 

*  OUTPUT  STARTING  VALUES 

WRITE  (6,490) 

WRITE  (6,500)  (PARAM(I) ,1=1,2) 

*  COMPUTE  IOGLIKELIHOOD  VALUE  FOR  INITIAL  STARTING  ESTIMATES 

CALL  LLIKE(l) 

FV=- L*FVAL 
WRITE (6, 510)  FV 

*  CHECK  FOR  CONVERGENCE 

CALL  CONCHK ( 1 ) 

IF  (CONFLG.EQ.l)  GO  TO  1000 

1  ITER=ITER+1 

WRITE  (6,520)  ITER 

*  IF  ITERATIONS  EQUAL  FIFTY,  PROGRAM  IS  TERMINATED  EARLY 

IF  (ITER. GE. 50)  THEN 
WRITE  (6,560) 

GO  TO  1000 
END  IF 

*  IF  IRES=1 ,  OBJECTIVE  FUNCTION  IS  NO  LONGER  IMPROVING  AND  GRADIENT 

*  IS  NOT  EQUAL  TO  ZERO.  THE  "A"  MATRIX  IS  RESET  TO  THE  IDENTITY  MATRIX. 

IF  ( ITER. EQ. IRES)  THEN 
A( 1 )=1 . 0D0 
A ( 2 )=0 . ODO 
A(3)=O.ODO 
A(4)=l .ODO 
IRES=IP.ES+20 
END  IF 

*  COMPUTE  DIRECTION  OF  IMPROVEMENT 

DIR( 1 )=-l*(A(i )*GRAD( 1 )+A(2 )*GRAD(2) ) 

DIR(2)=*1*(A(3)*GRAD(1 )+A(4)*GRAD(2) ) 

*  COMPUTE  LINE  STEP  SIZE  (ALPHA) 

CALL  LINE 

*  COMPUTE  DIFFERENCE  BETWEEN  PARAM (K+l )  AND  PARAM (K) ,  AND 

*  DIFFERENCE  BETWEEN  GRAD (K+l )  AND  GRAD(K) 

DO  10  1=1 ,2 

DELP ( I ) =PARAM ( 1+2 ) -PARAM ( I ) 

DELC-(l )  =GRAD  (1  +  2)  -  GRAD  ( I ) 
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CONTINUE 


*  UPDATE  'A'  MATRIX 


C1=DELP ( 1 ) *DELG ( 1 ) +DELP ( 2 ) *D2LG ( 2 ) 

A1 ( 1 ) =DELP ( 1 ) *DELP ( 1 ) /Cl 
A1  2)=DELP(l)*DELP(2  /Cl 
Al(3>=DELF(l)*DELP(2  /Cl 
A1(4)=DELP(2)*DELP(2)/C1 


C2=1-((DELP(1)*DELG(1) 
C3=-1*((DELP<1)*DELG(2 
C4=-l*( (DELP(2  *DELG(1 
C5=1-((DELP(2)*DELG(2) 


A2(1)=C2*(A(1)*C2+A(2)*C3)+C3*(A(3)*C2+A(4)*C3) 
A2(2  =C2*(A(1)*C4+A(2)*C5)+C3*(A(3  *C4+A(4)*C5 
A2(3)=C4*(A(1 )*C2+A  2)*C3)+C5*(A(3  *C2+A(4)*C3) 
A2(4)=C4*(A(1)*C4+A(2)*C5)+C5*(A(3)*C4+A(4)*CS) 

DO  20  1=1,4 

A(I )=A2( I )+Al (I ) 

20  CONTINUE 


*  UPDATE  OLD  PARAMETER  ESTIMATES  AND  GRADIENTS  TO  NEW  PARAMETER 

*  ESTIMATES  AND  GRADIENTS 


DO  25  1=1.2 

PARAM \ I )=PARAM ( 1+2 ) 

GRAD ( I ) =GRAD ( I + 2 ) 

25  CONTINUE 

*  OUTPUT  UPDATED  PARAMETERS 

WRITE  (6,500)  (PARAM(I) ,1=1,2) 

*  OUTPUT  FUNCTIONAL  VALUE 

FV=-1*FVAL 
WRITE (6 , 510)  FV 

*  CHECK  FOR  CONVERGENCE 

CALL  CONCHK(l) 

IF  (CONFLG.EQ . 1 )  GO  TO  1000 

*  START  NEW  ITERATION 

GO  TO  1 

*  PROGRAM  TERMINATION 
1000  WRITE  (6,530) 

WRITE  (6,540) 

WRITE  (6,550)  A 

*  OUTPUT  FORMAT 

490  FORMAT  ( 1 OSTARTING  VALUES : 1 ,/, 1  . ') 

£00  FORMAT  ( 1 OPHAT  IS:  ',F15.8,/,'  LHAT  IS:  '.F15.S) 

510  FORMAT  (' OLOGLIKELIHOQD  VALUE  IS:  ',F13.6) 

520  FORMAT  ( ' 0 ITERATION :  ',13,/,'  - - ' ) 

530  FORMAT  ( 1 OPROGRAM  COMPLETION') 

540  FORMAT  ( ' OFINAL  "A"  MATRIX  IS:1) 

550  FORMAT  ' 0 1 , 2 ( F15 . 8 ) , / , 1  ',2(F15.8)) 

500  FORMAT  (  '  OUNSUCCESSF’JL  TERMINATION:  ITERATIONS  =  50') 

STOP 

END 


SUBROUTINE  INFUT 


*  READS  FLIGHT  AND  MAINTENANCE  RECORDS  USED  IN  THE  CALCULATION 

*  OF  MIXED  MODEL  MAXIMUM  LIKELIHOOD  ESTIMATES  FOR  PHAT  AND 

*  LAMBDA  HAT 


IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/ RL/ ALPHA , FVAL ,UFVAL , DIR , FTNF , FTF , GRAD , PARAM 
COMMON/INT/CONFLG. COUNT, TFLAG, CENSOR, NF 
DIMENSION  CENSOR(200) ,NF(200) ,DIR(2) ,FTNF(200) 

DIMENSION  FTF (200) , GRAD (4) , PARAM(4) 

INTEGER  CENSOR, CONFLG, COUNT, TFLAG 

*  STARTING  POINTS  FOR  PARAM(l) ,PHAT,  PARAM ( 2 ) , LHAT ,  AND  RECORD  COUNTER 

PARAM(1)=.67D0 
PARAM(2)=.01D0 
COUNT* 1 

*  INPUT  FLIGHT  DA'iri 

5  READ  (3 , 500 ,END=10Q0)  CENSOR(COUNT) ,NF(COUNT) , FTNF (COUNT) , FTF(COUN 
CT) 

COUNT=COUNT+l 
GO  TO  5 

*  INPUT  FORMAT 

500  FORMAT  (II , IX, 14 , SX, F6 . 1 , IX, F3 . 1 ) 

1000  C0UNT=C0UNT-1 
RETURN 
END 


SUBROUTINE  CONCHK(N) 

*  DETERMINES  IF  OPTIMAL  SOLUTION  HA'S  BEEN  FOUND 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/ RL/ ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 
COMMON/INT/CONFLG , COUNT , TFLAG , CENSOR  ,NF 
DIMENSION  CENSOR ( ZOO) ,NF(2C0) ,DIR(2) , FTNF (200) 

DIMENSION  FTF ( 200 ), GRAD (4) , PARAM (4) 

INTEGER  CENSOR. CONFLG, COUNT, TFLAG 
REALMS  NORMG 

*  TERMINATE  IF  THE  NORM  OF  THE  GRADIENT  (NORMG)  IS  LESS  THAN  THE 

*  THE  SPECIFIED  CRITICAL  VALUE  (CRIT) 


NORMG* ( ( GRAD ( N ) *GRAD (N) )  +  ( GRAD ( N+ 1 ) ’’'GRAD (N+l  )))**. 5 

IF  (NORMG. LE.l.)  THEN 
CONFLG=l 
WRITE  (6,500) 

END  IF 

FTEST=ABS (UFVAL-FVAl ) 

IT  (FTEST.LE. . 000001. AND. NORMG. GT.l.)  THEN 
CONFLG* 1 
WRITE  (6,510) 

END  IF 

*  OUTPUT  FORMAT 

500  FORMAT  ( 1 ONORMAL  TERMINATION  -  OPTIMAL  FOUND1) 

510  FORMAT  ( 1 0 PROGRAM  TERMINATION  -  LACK  OF  FUNCTIONAL  IMPROVEMENT') 

RETURN  , 

END  i 


I 
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SUBROUTINE  LINE 


COMPUTES  MAXIMUM  STEP  SIZE  (ALPHA) 


IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/ RL/ ALPHA , FVAL , UFVAL ,  DIR , FTNF , FT F , GRAD , PARAM 
COMMON/  INT/CONFLG ,  COUNT ,  TFLAG ,  CENSOR ,  NF 
DIMENSION  CENSOR (200) , NF(200 ) ,DIR(2), FTNF ( 200) 

DIMENSION  FTF (200 ) , GRAD (4) , PARAM(4) 

INTEGER  CENSOR ,  Cf'NFLG ,  COUNT ,  TFLAG 
REAL*8  LB 

UFVAL-FVAL 

K=1 

TFLAG=0 
ALPHA=O.ODO 
ALPHAK=.0001DO 
LB=0 . ODO 

COMPUTATION  OF  SEARCH  INTERVAL  TAU.  TAU  IS  APPROXIMATELY 
.618  AND  TAU1  IS  APPROXIMATELY  .382 
TAU=(5**.5-l)/2 
TAU1=1-TAU 

BOUNDING  PHASE  -  SWANN'S  METHOD 
STEP- . OOOOIDO 

PARAM(3 )=PARAM( 1 ) -5TEP*DIR( 1 ) 
PARAH(4)=PARAM(2)-STEP>'DIR(2) 

CALL  LLIKE(3) 

F1=FVAL 

P AhAM ( 3 ) =PARAM ( 1 ) +STEP*DIR ( 1 ) 
PARAM(4)=PAPAM(2)+STEP*DIR(2) 

CALL  LLIKE ( 3  )  • 

F2=FVAL 

IF  (F1.GE.F2)  THEN 
F1=F2 
GO  TO  5 
END  IF 

IF  (F2.GT.F1)  STEP=~1*STEP 

ALPHA=ALPHA+ ( 2**K ) *STEP 
PARAM ( 3 ) =P ARAM ( 1 ) +ALPHA*DI R ( 1 ) 

PARAM ( 4 ) =PARAM ( 2 ) +ALPHA*D IR ( 2 ) 

COMPUTE  LOGLIKELIHOOD  VALUE 
CALL  LLIKE (3) 

F2=FVAL 

IF  NEW  VALUE  IS  LESS  THAN  PREVIOUS  VALUE  TERMINATE  BOUNDING 

IF  (F2.GT.F1)  GO  TO  10 

F1=F2 

K=K+1 

GO  TO  5 

UB=ALPHA 

CONDUCT  GOLDEN  SECTION  SEARCH 

B1=TAU1*(UB-LB)+LB 

PARAM 3  )=PARAM(  1  )+El*DIR(l) 

PARAM(4)=?ARAM(2)+B1*DIR(2) 

CALL  LLIKE (3) 

FF 1  =CVAT 

B2=TAU*(UB-LB)+LB 

PARAM (3  -PARAM (1 )+B2*DIR(I) 

PARAM ( 4 )=PARAM ( 2 ) +32>DIR ( 2 ) 

CALL  LLIKE ( 3 ) 


*  * 


50 


FB2»FVAL 

IF  (FBI .LT.FB2)  THEN 
ALPHA=B1 
UB=B2 
B2=B1 
FB2*FB1 

B1«TAU1*(UB-LB)+LB 
PARAM(3)=PARAM(l)+Bl*DIR(l) 
PARAM(4)=PARAM(2)+B1*DIR(2> 
CALL  LLIKE(3) 

FB1=FVAL 
GO  TO  100 
END  IF 

IF  (FB2.LE.FB1)  THEN 
ALPHA*B2 
LB=B1 
B1=B2 
FB1=FB2 

B2=TAU* ( UB-LB ) +LB 
PARAM(3  =PARAM(1)+B2*DIR(1) 
PARAM ( 4 ) =PARAM ( 2 ) +B2*DIR( 2 ) 
CALL  LLIKE (3) 

FB2=FVAL 
GO  TO  100 
END- IF 

*  CHECK  FOR  LINE  SEARCH  TERMINATION 
100  CALL  TERM(ALPHAK) 

ALPHAK=ALPHA 

IF  (TFLAG.EQ.l)  GO  TO  5000 
GO  TO  50 

5000  PARAM  ( 3 )  =PARAM(  1 )  +ALPHA'tDIR  ( 1 ) 
PARAM(4)=PARAM(2 )+ALPHA*DIR( 2  > 
CALL  LLIKE(3) 

RETURN 

END 


SUBROUTINE  TERM(ALPHAK) 

*  DETERMINES  IF  CONVERGENCE  CRITERIA  IS  ESTABLISHED  FOR  LINE 

*  SEARCH  TERMINATION 


IMPLICIT  REAL*8  (A-Hr0-Z) 

COMMON/RL/ ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 
COMMON/ INT/CONFLG , COUNT , TFLAG , CENSOR , NF 
DIMENSION  CENSOR( 200 ) ,NF(2Q0) , DIR ( 2 ) , FTNF( 200 ) 
DIMENSION  FTF (200) , GRAD (4) ,PARAH(4) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

TFLAG=0 

CTEST=ABS ( (ALPHAK- ALPHA) /ALPHA) 

IF  (CTEST.LE. .003)  TFLAG=1 

RETURN 

END 


SUBROUTINE  LLIKE(N) 

COMPUTES  LOGLIKELIHOOD  FUNCTIONAL  VALUE  AND  GRADIENTS  FOR 
FOR  THE  WEIBULL  MODEL 

IMPLICIT  REAL*8  (A-H.O-Z) 
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INPUT  FLIGHT  DATA 
DO  5  1=1 , COUNT 


IF  (FTF(: 
IF  (FTNF 
Cl*  ( FTNF 
C2*PARAM 
C3=PARAM 


FTF(I).E0.0.)  FTF(I)=. 000001 
FTNF ( I ) . SO . 0 . )  FTNF(I)*. 000001 


)  .EO.O.) 
+FTF(I)) 
*C1 

)*FTNF(I) 


C4=EXP ( - 1* (C2**PARAh ( 1 ) )) /EX? ( - 1*( C3*+?K 
SUMKI1=SUMKI1+C* ( ( ( (FTNF( I ) **PARAM( 1 ) -Cl 

*PARAM( 1 ) -C1**PARAM( 1 ) )*C4) )/ ( l-C4)+< 
)**PARAM(l))*EXP(-l*?C3’,«*PARAM(l)n) 
SUMKI2=SUMKI2+(  ('( (C3*PARAM(1 ) *L0G(C3)) -  ( 


-1*<C3**PARAM(1))) 

ARAM ( 1 ) -Cln*PARAM( 1 ) )*C4 


SUMKI2=SUMKI2 
C  **2)/{l-C4) 
CONTINUE 


C3**PARAM(l)*LOG(C3) ) 


Cl  **PARAM(  1 ) ) *C4 ) * ( < FTNF  ( I )  * 
+ (FTNF ( I )**PARAM( 1) )*  ( FTNF ( I 

- (C2**PARAM( 1 )*LOG(C2) ) )*C4) 
2 


PSD1* 

PSD2= 


*SQRT ( 1/SUMKI1) 
;=SQRT  ( 1/SUMKI2) 


*  CONFIDENCE  INTERVALS 

LL1=PARAM( 1 )  -1 . 96*PSD1 
UL1=PARAM(1 )*1 . 96*PSDi 
LL2*PARAM(2)-1.96*PSD1 
UL2=PARAM(2)t-l.96’‘TSDl 

*  OUTPUT 

WRITE  (6,100)  PSD1 ,LL1 ,UL1 
WRITE  6,110  PSD2 ,LL2 ,UL2 
100  FORMAT  ( ' OPSDl i  '^6.5,'  CONFIDENCE  INTERVAL 

110  FORMAT  ( 1 0PSD2 :  ' , F6 . 5 , 1  CONFIDENCE  INTERVAL 

RETURN 

END 


(',F6.5,',',F6.5,'V) 

(',F6.5,',',F6.5,')') 


SUBROUTINE  LLIKE(N) 

*  COMPUTES  LOGLIKELIHOOD  FUNCTIONAL  VALUE  AND  GRADIENTS  FOR 

*  THE  MIXED  MODEL. 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON /RL/ ALPHA ,  FVAL  ,  UFVAL  ,  DIR,  FTNF ,  FTF  ,GRAD  , PARAM 
COMMON/ INI/ CONFIG, COUNT, TFLAG, CENSOR, NF 
DIMENSION  CENSOR(200) ,NF(200) ,DIR(2) ,FTNF(200) 

DIMENSION  FTF(200) ,GRAD (4) , PARAM (4 ) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

*  ESTABLISH  BOUNDARY  CONDITIONS 

IF  (PARAM(N) .LT. .0000000001)  PARAM (N) =. 0001 
IF  PARAM  N.GE.l.)  PARAM (N)  =  . 9999 
IF  (PARAM(N+1) .LT. .0000000001)  PARAM(N+1 )= . 0001 

FVAL=0 . 

GRAD(N)=0 . 

GRAD(N+1)=0. 

Cl=l -PARAM (N) 

C2=L0G(C1) 

*  INPUT  FLIGHT  DATA 

DO  5  1=1 , COUNT 

*  UNCENSORED  RECORDS 

IF  (CEMSOR(I) .EQ.O)  THEN 

C3=EXP(*-PARAM(N-M )  *FTF(  I ) ) 

*  FUNCTIONAL  VALUE  CALCULATIONS  ^  „ ,  T . 

FVAL=FVAL-LOG( 1-C3*C1 ) - ( NF ( I ) - 1 ) *C2+PARAM (N+l ) *FTNF ( I ) 

*  GRADIENT  CALCULATIONS  ^  ,  ,  ,  v  w 

GRAD(N)=GRAD(N)-(C3/(1-C1*C3))+(NF(I)-1)/C1 

GRAD(N+1)*GRAD(N+1)-(FTF(I)*C1*C3)/(1-C1*C3)+FTNF(1) 


»  * 


I  ;  »-  -’.««>  .'-fcs  «',<  '»*,*  *«>•.«■»«  ».«*.<  *-*  •,»  *«i.i  y  *.«  »,  Lr*v  i  < 


GO  TO  5 
END  IF 

CENSORED  RECORDS 

FUNCTIONAL  VALUE  CALCULATIONS 

ORU,I^LSLCuSUo;iSC2*!'A!i"<Nn  >  *™F(  1  > 

GRAD(N)=GRAD(N)+NF(T.)/C1 
GRAD(N+1 ) =GRAD (N+l ) +FTNF ( I ) 

CONTINUE 

RETURN 

END 
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APPENDIX  H 

WEIBULL  SIMULATION  OF  COMPONENT  FAILURES 
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VARIABLE  DEFINITIONS 
ALPHA.  WEIBULL  SHAPE  PARAMETER 

AC. MATRIX  USED  TO  STORE  EACH  AIRCRAFTS  PREVIOUS  FLIGHT  HOURS 
WITHOUT  FAILURE  AND  UNCONDITIONAL  SURVIVAL  FUNCTION 
AS .AIRCRAFT  SELECTED  FOR  FLIGHT  EVOLUTION 
CPF i CONDITIONAL  PROBABILITY  OF  FAILURE 
DTIME: DEPLOYMENT  TIME 

FAIL .MATRIX  USED  TO  STORE  SIMULATED  FAILURES  AT  SPECIFIED  TIME 
INTERVALS  FOR  EACH  DEPLOYMENT  SIMULATION 
FT .FLIGHT  TIME 

GFC .GENERATED  FAILURE  CRITERIA 

HM .HISTORICAL  MEAN 

HSD. HISTORICAL  STANDARD  DEVIATION 

LHAT .MAXIMUM  LIKELIHOOD  ESTIMATE  FOR  WEIBULL  SCALE  PARAMETER 
NAC .NUMBER  OF  AIRCRAFT  USED  IN  SIMULATION 

NI .NUMBER  OF  TIME  INTERVALS  TO  BE  USED  IN  SIMULATION  (500  HR  INT.) 
NS .NUMBER  OF  SIMULATIONS  TO  BE  RUN 

NSEED1  .SEED  USED  TO  DETERMINE  AIRCRAFT  CHOSEN  FOR  FLIGHT 
NSEED2  .SEED  USED  TO  DETERMINE  FLIGHT  DURATION 

NSEED3  .SEED  USED  TO  HDETERMINE  IF  A  COMPONENT  FAILURE  WAS  OBSERVED 

TGATE  .'TIME  GATE 

TIME .ACCUMULATED  FLIGHT  TIME 

UCPS  .UNCONDITIONAL  PROBAIMTY  OF  FAILURE 

DIMENSION  AC (70 ,  2) ,  FAIL(20 , 100) , EST ( 20 , 4 ) 

INTEGER  AS 
REAL  *  4  LHAT 

INITIALIZE 

DATA  LHAT/. 0149 /ALPHA/. 807 9/ 

DATA  NI/16/NS/100/NAC/48/HM/2. 18/HSD/ .9897/ 

DATA  NSEED1/749684/NSEED2/4683957/NSEED3/ 5872654/ 

COMMENCE  SIMULATION ( I ) 

DO  1  1=1, NS 
NF=0 
TIME=0. 

TGATE=5GQ . 

DO  5  K=1 ,NAC 
AC (K, 1 ) =0 . 

AC(K,2)=1. 

CONTINUE 

SIMULATE  AND  RECORD  COMPONENT  FAILURES  USING  500  HR  INTERVALS 
DO  10  J=1 ,NI 

LRND  AND  LNORM  ARE  SUBROUTINES  IN  A  RANDOM  NUMBER  GENERATOR  PACKAGE 
CALLED  LLRANDQMI I .  A  NON-IMSL  PROGRAM  AVAILABLE  AT  THE  NAVAL 
POSTGRADUATE  SCHOOL. 


GENERATE  UNIFORM  (0,1)  RANDOM  VARIABLE  TO  DETERMINE  WHICH  AIRCRAFT 
WILL  BE  FLOWN 

CALL  LRND (NSEEDl.U, 1,2,0) 


AS=(U*(NAC-1 ) }+l . 5 

GENERATE  NORMAL  (0,1)  RANDOM  VARIABLE  TO  COMPUTE  FLIGHT  DURATION 


*  * 


gun wawa W HBB BU law  W  TOWyJWW rv*  r« r^nj^-A r>»a*-aji inan  - j  .  J  j™r/  ruvl'u.»’wi«  ^« 


CALL  LN0RM(NSEED2 ,RN,  1,2,0) 

FT® ( RN*HSD ) +HM 
IF  (FT.LT. .5)  FT=  •  5 
IF  (FT.GT.6.)  FT=6. 

*  COMPUTE  THE  UNCONDITIONAL  PROBABILITY  OF  SURVIVAL  FOR  S+FT  HRS, 

*  WHERE  S  EQUALS  THE  TOTAL  FLIGHT  HOURS  FLOWN  SINCE  THE  LAST 

*  COMPONENT  FAILURE,  AND  FT  EQUALS  THE  FLIGHT  TIME  OF  THE  SORTIE 

*  BEING  FLOWN. 

UCPS=EXP( -1*( (LHAT*(AC(AS , 1 )+FT) )**ALPHA) ) 
CPF=1-UCPS/AC(AS,2) 

GENERATE  UNIFORM  (0,1)  RANDOM  VARIABLE  TO  DETERMINE  IF  COMPONENT 
FAILURE  OCCURS 

CALL  LRND(NSEED3,GFC, 1,2,0) 

*  COMPONENT  FAILURE 

IF  (GFC.LT.CPF)  THEN 
NF=NF+1 
AC( AS , 1 ) =C . 

AC(AS,2)=1. 

GO  TO  15 
END  IF 

*  COMPONENT  SURVIVES 

AC(AS  ,  1 )=AC( AS , 1 )+FT 
AC(AS , 2)=UCPS 

*  UPDATE  TIME 

15  TIME=TIME+FT 

IF  (TIME.GE.TGATE)  THEN 
FAIL( J , I )=NF 
TGATE=TGATE+500 . 

GO  TO  10 
END  IF 

*  NEXT  AIRCRAFT  SORTIE  TO  BE  SIMULATED 

GO  TO  11 

*  NEXT  TIME  INTERVAL 

10  CONTINUE 

*  NEXT  DEPLOYMENT 

1  CONTINUE 

*  CALCULATE  MEAN  AND  STANDARD  DEVIATION  FOR  THE  NUMBER  OF  FAILURES 

*  OBSERVED  AT  EACH  TIME  INTERVAL  (MEAN=EST( J, 1 ) ; STANDARD 

*  DEVIATION=EST ( J , 2 ) ) 

DO  20  J=1,NI 
DO  25  1=1, NS 

EST(J,1)=EST(J,1 )+FAIL(J,I) 

F^T  f  7 , 2 )=EST' 7,2; TAILVJ  T  >*FAIL( J , I ) 

2'  CONTINUE 

20  CONTINUE 


1!t( (E^t{ (N5*EST( J , ] )*EST(J , 1) ) >/ (NS-1) )**. 5 
30  CONTINUE 

*  95  PERCENT  CONFIDENCE  INTERVAL  (LOWER  LIMIT=EST( J , 3) ; 

*  UPPER  LIMIT=EST( J,4) 

D)  35  J=1 ,NI 

HINT=1.96*EST(J,2) 

EST(J,3)-EST(J,1)-HINT 
EST(J,4)=EST(J,1)+HINT 
35  CONTINUE 

*  OUTPUT 

TIME=500. 

*  WRITE  (6,1000) 

DO  40  J=1  NI 

WRITE"  (6,1100)  TIME ,EST(J,1),EST(J,2) ,EST( J ,3) ,EST( J ,4) 
TIME=TIME+SOO. 

40  CONTINUE 

*  OUTPUT  FORMAT 

*1000  FORMAT  ('1  TIME  MEAN  STD  DEV  LOWER  UPPER') 

1100  FORMAT  (F5.0,4(2X,F7 .3) ) 

STOP 

END 
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GRAPHICAL  COMPARISONS 
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Figure  1.5  Flight  Hour  Models:  Model  Set  WLCMOX 1600. 
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Figure  1.6  Geometric  Model:  Model  Set  \VLC:-U>X  |(>W. 
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Figure  1.10  Geometric  Model:  Model  Set  W  IX':4622100. 
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Figure  1.12  Geometric  Model:  Validation  Set  \VL’C:4622IOO. 
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Figure  1.16  Geometric  Model:  Validation  Set  \VL'C:5o\250U. 
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Figure  1.18  Geometric  Model:  Model  Set  WUC:(02Zlou 


Figure  1.19  Flight  Hour  Models:  Validation  Set  WLC:632Z1<>U. 
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Figure  1.22  Geometric  Model:  Model  Set  WL'CMA  1500. 
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Figure  1.24  Geometric  ?vl odcl:  Validation  Set  \\'l'C:74AI5Ui>. 
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Figure  1.26  Flight  Hour  Models:  Validation  Set  WL'C:74A5\lW. 
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Figure  1.31  Geometric  Model:  Validation  Set  \V U C : 5 7 7 2 
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Figure  1.32  Flight  Hour  Models:  Model  Set  \YLC:091SIW 


Figure  1.34  Flight  Hour  Models:  Validation  Set  WL'C.^IS  100 
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